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Preface 


This technical report is an expansion of “Higher order asymptotic corrections applied in an 
EM algorithm for estimating educational proficiencies”. It includes a brief overview of the 
MGROUP computer program originally written by Robert Mislevy and Kathy Sheehan. 
The primary purpose of this report is to provide detailed documentation of the E-step 
calculations implemented for the 1992 NAEP assessments. This report does not describe 
the input files and commands required to execute the MGROUP program. 


Abstract 


Mislevy(1984, 1985) introduced an EM algorithm for estimating the parameters of a 
latent distribution model that is used extensively by the National Assessment of Educational 
Progress. Second order asymptotic corrections are derived and applied along with more 
common first order asymptotic corrections to approximate tke expectations required by the 
E-step of the EM algorithm. These corrections produce much more accurate approximations 
than those produced by two different normal approximations. The asymptotic corrections 
are applicable in high dimensional models. The number of function evaluations required by 
these methods grows linearly as dimension iucreases, in contrast to the exponential growth 


with dimension of product quadrature rules which yield similar accuracy. 


Keywords and phrases. EM algorithm, higher order asymptotic corrections, Laplace 


expansion, item response models. 
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1 Introduction 


A multivariate multiple regression model with p latent variables representing proficiency in 
several closely related subject areas is used extensively to summarize educational testing 
data obtained from national surveys. For a sample of n examinees, the proficiency variables, 
denoted by 6; = (6;1,..., Op), 1 = 1,...,n, are assumed to follow a multivariate normal dis- 


tribution, 0; ~ N (I'z;, Z), where z = (zi1,..., Zim), X = (21,...,2n)', are m fully mea- 


sured demographic and educational characteristics for each examinee, I = [y ere | | 


are the m*p parameters of the regression of 8 on x, and Z is the variance-covariance matrix 
for @ conditional on x. The (0;,...,4,) are often referred to as scales in the psychomet- 
ric literature. For example, several large mathematics assessments use the scales numbers 
and operations, measurement, geometry, and algebra. This model is used extensively for 
scaling, equating, generation of multiple imputations, and other reporting purposes by the 
National Assessment of Educational Progress (NAEP), (Johnson and Zwick, 1990, Mislevy, 
1991, Mislevy, Johnson, and Muraki, 1992), and the Adult Literacy Survey, (Kirsch and 
Jungeblut, 1986). An overview of he statistical methods used in NAEP is given in Beaton 
and Zwick(1992). 

The variable 0; is not measured directly, but information about @; is obtained from 
several examination items given to each examinee denoted by y;, Y = (y;,...,Y,). Con- 
ditional on the 0;, the responses of different examinees are assumed to be independent 
with a distribution denoted by []/L, f (y; | 9:). The key feature used to form improved 
asymptotic approximations is that each examination item contributes information about a 
single scale. Decomposing the responses to the examination items for the i" examinee as 
¥;= (var, uP ip): corresponding to the scales (6;1,...,6:p) respectively, f (y,; | @;) can be 


written in the form 
f (ys 181) = fa (win | 811) => Sy (Yip | Op) = (81) = bp (8p) S18). (1.1) 


Conditional on the @,;, the responses to the examination items measuring the same 


proficiency scale of an examince are also assumed to be independent, so that f, (v.,) is 


] 


the product of the response probabilities tor each examination item administered. For 
dichotomous items represented by y = 0,1, a three parameter logistic item response model 
is used, 

P(y=1|6)=7+(1-—7)/{1+exp[a(@-6)]} , (1.2) 
and for ordered polytomously scored items, y = 1,...,c, a generalized partial credit model 


is used, 
exp [Dj-1 (9 - 4;)| 


Pa exp (het a(6- 6;)| 
where §,; = 0. Hambleton and Swaminathan(1985), and Bock and Aitkin(1981) describe 


P(y=k|@)= (1.3) 


the model and the estimation of the item specific parameters (a, 6,7) for the dichotomous 
responses, and Muraki(1992) describes the polytomous response model. 
With the item specific parameters estimated separately and regarded as known, the 


likelihood function for (I", Z') based on the (X,Y) data is given by 


TT] [fi (ui | 81) +> fo (vip | in) M(Ois T’m, E) a8; (1.4) 
t=1 


The function in (1.1) is called the examinee likelihood function because it is the likelihood 
function for estimating 0; from the y; data, and the integrand in (1.4) is called the examinee 
posterior distribution (conditional on (I, X) ). 

Mislevy(1984, 1985) formed an EM algorithm, Dempster, Laird, and Rubin(1977), to 
estimate (I, ) by considering the 0; as missing. If 0; were observed, the maximum 


likelihood estimators would be 


6; 

P= [4.14] = (AX) Xt |, (1.5) 
6, 

s-= s=15 (6; - i'z;) (8; - f'z,) (1.6) 


t=] 


When the sufficient statistics for the complete data estimators in (1.5) and (1.6) are re- 
placed by their expected values conditional on the data and provisional parameter values, 


(I, £&;), straightforward calculations similar to those in Little and Rubin(1987) show that 


the updated parameter estimates are 


Pin = (XX) X'| : |, (1.7) 
a, 
1{2 et - ; 

Dun = = [power IXY. PE) +> (8; - Pisizi) (6 - Pis121) - (18) 
where 6; = (Bi1,... 5p), = E(0; | X,Y,I;, Z:) is the examinee posterior mean. This 
process is repeated until convergence. 

Equation (1.8) decomposes 41 into a component representing the remaining uncer- 
tainty in 0; after observing (z;,y,;) and the variation of the estimates of 0; about the 
regression of @ on X. When the examination items, y,;, determine 6; so that 6; = 0;, there 
is no remaining uncertainty about 0; and the EM algorithm becomes a standard multivari- 
ate multiple regression. In the other limiting case when y; contains no information about 
6;, then 0; = F'izi, Pia) = Fy, and Dy4) = Zi. 

An outline of the paper follows. Section 2 presents an example using NAEP data with 
emphasis on the difference in the results obtained using different numerical approximations 
for the expectations in the E-step of the EM algorithm. Section 3 describes a method for 
computing the examinee posterior moments required for the E-step using a common first 
order asymptotic correction for the examinee posterior means, and a specialized second 
order asvmptotic correction derived for the examinee posterior variance-covariance matrix. 
Section 4 describes the two alternative approaches illustrated in Section 2 for computing 
the examinee posterior moments. Section 5 summarizes the relative computing speed of the 
different methods. Section 6 contains a brief overiveiw of the MGROUP program. Sections 


7, 8, and 9 give a detailed description of the different E-step algorithms. 
2 Example 


Table 1 presents results using data from the main NAEP sample of 8790 students for the 
1990 assessment of mathematics for students in the fourth grade or nine years of age. Results 


are presented for the Numbers and Operations scale, in which examinees are given between 


3 


19 and 26 examination items, and the Measurement scale, in which examinees are given 


between six and eleven examination items. An example with only two scales is presented 
because it is possible to apply nearly exact numerical quadrature in this setting. The 
example uses a subset of the complete set of NAEP conditioning “X” variables. There are 
eleven dichotomous variables representing additive factors for sex, race, modal age-grade 
status, and teacher rating of class ability. A description of the covariate labels in Table 1 
is given in Table 2. 

The sample size for typical NAEP applications of the model in Section 1 in the various 
surveys range from 1000 to 10,000 examinees, with one to six scales. The EM algorithm 
often requires 100 or more iterations to reach convergence, thus requiring the evaluation of as 
many as (10,000 examinees) * (100 iterations) * (6 means + 21 variances and covariances) 
= 27 million six dimensional integrals. As a consequence, methods for evaluating these 
integrals must be very efficient. 

Table 1 displays the converged estimates of (I", 3’) from the EM algorithm based on 
three different methods of evaluating the examinee posterior moments. The estimates in 
the first two columns of Table 1 are based on a simple numerical quadrature method for 
the E-step described in Sections 4 and 7. By using additional quadrature points, and 
simulatiou >tudies, the numerical accuracy of this method has been confirmed. However, 
the quadrature approach is very slow in applications with two scales (four to eight hours for 
a typical apnlication on a Sun SPARC Station IPX), and computationally unfeasible when 
there are more than two scales. 

The second set of estimates displayed in the middle columns of Table 1 are based on an 
Integrated normal approximation motivated by the asymptotic normality of the examinee 
likelihood function, /(@), as the number of examination items increases. This approach, 
which is described in Sections 4 and 8, produces a very fast algorithm useful for checking 
purposes and obtaining starting values. For the number of examination items commonly 
administered, however, the accuracy of this algorithm is inadequate. Discrepancies between 


the nearly correct estimates obtained using numerical quadrature and the estimates ob- 


i0 


Table 1: Parameter estimates for the mathematics example 


Numerical Integrated Normal Higher Order 
Quadrature Approximation Asymptotic Corrections 


Num&Op Measurement Num&Op Measurement | Num&Op Measurement 


Intercept -0.569 -0.357 -0.561 (0.03)  -0.362 (0.04) -0.568 -0.354 


Female -0.026 0.114 -0.021 (0.02) -0.138 (0.02) | -0.026 0.115 
Afric -0.512 0.643 -0.552 (0.02) -0.828 (0.03) | -0.516 0.651 
Hisp -0.407 0.442 -0.441 (0.02) -0.566 (0.03) | -0.410 0.448 
Asian 0.075 0.041 0.070 (0.05) 0.038 (0.06) | 0.074 0.042 
EMALMG | -0.577 0.547 -0.614 (0.03) -0.700 (0.04) | -0.582 0.552 
EMAGMG | __ 0.572 0.483 ° | 0.557 (0.16) 0.522 (0.21) | 0.573 0.484 
GMAEMG | -0.193 0.156 -0.200 (0.02) -0.208 (0.03) | -0.194 - 0.157 
LMAEMG | 0.140 -0.356 0.155 (0.15)  -0.338 (0.20) | 0.143 0.358 
CLABH 0.441 0.407 0.447 (0.04) 0.467 (0.06) | 0.443 0.406 
CLABM 0.012 0.014 0.019 (0.03) 0.022 (0.03) | 0.012 0.013 
CLABL -0.190 -0.223 -0.203 (0.04)  -0.273 (0.05) | -0.190 0.224 

0.364 0.344 0.382 0.456 0.368 0.350 
z= 0.344 0.368 0.456 0.563 0.350 0.374 


Note: The numbers in parentheses are standard errors (SE). These SE are also used 
as approximate SE for the other computational methods. 


Table 2: Description of regressor variables 


Regressor | Description 

Intercept 

Female One if examinee is female 

Afric One if examinee is African-American 

Hisp One if examinee is Hispanic 

EMALMG | Equal to modal age, less than modal grade 

EMAGMG | Equal to modal age, greater than modal grade 
‘ GMAEMG | Greater than modal age, equal modal grade 

LMAEMG | Less than modal age, equal modal grade 

CLABH Teacher rates class ability high 

CLABM Teacher rates class ability medium 

CLABL Teacher rates class ability low 


cr 


tained using the normal approximation in Table 1 are large in substantive terms and when 
compared to approximate standard errors of the estimates. These approximation errors do 
not decrease as the number of examinees in the sample increase. The approximation errors 
are larger for the mecsurement scale because it is estimated from fewer examination items 
per examinee. The estimate of the residual variance for this scale is especially poor. The 
magnitude of these approximation errors are typical of applications with only a few exam- 
ination items per examinee. The errors in the parameter estimates are a consequence of 
errors in the approximation of both the examinee posterior means and variances in equations 
(1.7) and (1.8). The more common modal normal approximation to the examinee posterior 
distribution, described in Section 3, produces worse approximations than the integrated 
normal approximation so results for it are not presented in Table 1. 

The estimates presented in the final columns of Table 1 are based on the E-step examinee 
posterior moments computed using Laplace expansions, (Kass and Steffey, 1989, Mosteller 
and Wallace, 1964, and de Bruijn, 1981). The calculation of first order corrections for the 
6;, and the calculation of second order corrections for var (0; | X, ¥,I7:, Zt) are relatively 
simple and very efficient for the special structure of the examinee posterior distributions 
given in (1.1) and (1.4). The results in Table 1 are indicative of the good accuracy obtained 
from this approach even in difficult settings where examinees receive between five and ten 
items per scale. Details of this approach to evaluating the E-step integrals are given in 
Sections 3 and 9. 

The asymptotic corrections achieved the marked improvement by producing much more 
accurate approximations of the integrais computed for each examinee. To demonstrate 
this fact, the posterior means and variances were computed for each examinee using the 
three different methods with the same values of I, and ZX equal to the converged value 
of I and Z based on the numerical quadrature method. The root mean squares of the 
differences between the means, variances, and the covariance computed for each examinee 
using the asymptotic corrections and the numerical] quadrature range from 90 to 95 percent 


smaller than the corresponding differences between the Integrated normal approximation 


and the numerical quadrature. The performence of the Integrated normal approximation 
was especially poor for the variance of the shorter measurement scale. Figure 1 displays 
a plot of the posterior variances of 500 randomly selected examinees computed using the 
Integrated normal approximation versus the posterior variances computed using numerical 
quadrature. The Integrated normal approximation over-estimates the larger variances. This 
approximation error inflates the estimates of Z and it also produces a tendency for the 
magnitude of the I’ estimates to be too large as the EM algorithm iterates because the 
subsequent examinee posterior means are not shrunk sufficiently. Figure 2 shows that 
the asymptotic corrections dramatically reduce the apprcximation errors in the examinee 


posterio: variances. 


3 Higher order asymptotic corrections for the examinee 
posterior moments 


The EM algorithm in (1.7) and (1.8) requires the computation of the examinee posterior 
mean 6 = E(@| X,¥,I1, £1) and var(@ | X,¥,I:, 2). A computationally efficient ap- 
proach to approximating these moments uses asymptotic approximations as the number of 
examination items increases. Denote the negative log posterior distribution for an examinee 
with covariates, x, and examination item responses, y, by h (0) = — log|f (y | 0) N(0; Iz, 
~;)], and the total number of items by K. For a sequence of examination items increasing 
proportionately for each scale, the common modal (Bayesian) asymptotic approximation for 
the posterior mean and the variance-covariance matrix based on Laplace’s method, Kass 


and Steffey(1989), has the form 


6, = 6;+O(1/K) , (3.1) 
cov (6), 6% | X,Y L141, Zt41) = Gye + O(1/K?) , (3.2) 


a -1 a 
where 6; is the posterior mode of 6; and Gjx = ( weba: ase . The 6, are found by Newton- 
Raphson maximization requiring between two and three iterations. A rigorous statement 


and derivation of (3.1) and (3.2) and the asymptotic claims that follow are more difficult. 


for educational measurement models because examinees do not have identically distributed 


responses to differing examination items. Modifications of the usual asymptotic derivations 
that account for the varying characteristics of the examination items and the estimation of 
item specific parameters, like those in Haberman(1977), are not developed here. 

Higher order asymptotic corrections are obtained by computing the leading omitted 
terms in (3.1) and (3.2). For the special structure of h implied by (1.1) and (1.4), the 


computation of the correction terms is greatly simplified because the examinee likelihood 


function in (1.1) factors into distinct functions of (4),...,0,) so that 
Olog [1 (@)] _ 
00;00, ds 


and likewise for all higher order mixed derivatives. In addition, the normal prior distribution 
for 8 implies that all third and higher derivatives of the log prior distribution are zero, so that 
the only non-zero derivatives of third and higher order of h are the pure partial derivatives 
denoted by Ne j=1,...p, and r = 3,4,..., where the hat indicates that the derivatives 
of h are evaluated at 8. 
Using the special form of the derivatives of h, the first order correction to the mean in 
(3.1) is given by 
0; = 6, - 5 SG ieGerh, (3.3) 


r=] 


which is a specialization of the first order correction formula given by several authors and 
summarized in Kass and Steffey(1989). The higher derivatives of the log of the examinee 
likelihood function for the item response models in (1.2) and (1.3) can be computed effi- 
ciently, and these derivatives can be reused for the moments of each scale. Because only pure 
partial derivatives are required in (3.3), a direct approach to computing asymptotic correc- 
tions is more efficient in this setting than alternative approaches that avoid the calculation 
of higher derivatives such as Tierney, Kass, and Kadane(1989). 

To improve the usua] asymptotic variance-covariance matrix approximations in (3.2), the 
second order correction term was derived using the symbolic algebra program Mathematica. 


The Laplace expansion was applied to the integrals represented in standard form, Tierney. 
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Kass, and Kadane(1989) yielding the approximation 


cov as | X, Y, P41, £41) 
= Gie- DS (GirGerGrr) AM + = sEb| E - iG = 5) APADG,g 


21 s=1 


x (GjsGesGrr + GjsGkrGrs + GjrGksGrs + GjrGkrGss) } . (3.4) 


The derivation of (3.4) for the case with only one scale is presented in Appendix A, and 
extensions to the multi-scale case are indicated. 


When there is only one scale, (3.4) specializes to 
1 . o 2 
ar (6; | X,¥, Pes, Desa) ¥ Gu — sGh hi + Gt (A?) (3.5) 


As anticipated, the correction does not depend on the sign of A®), and it increases the vari- 
ance when the examinee posterior density decays slower than the normal distribution in the 
neighborhood of 6. When applied to the beta, gamma, and t, families of distributions, (3.5) 
produces corrections exactly equal to the corrections based on the expansion of the moment 
generating function in Tierney, Kass, and Kadane(1989). A simple extension of their The- 
orem 3 shows that (3.5) and the moment generating function approach produces identical 


variance formulas in one dimensional problems whenever the expansions are applicable. 


4 Alternative methods of computing the examinee poste- 
rior moments 


Numerical quadrature provides an accurate and computationally feasible approach for eval- 
uating the examinee posterior moments required by the E-step for applications with one 
or two scales. The quadrature implemented for applications with one scale evaluates the 
examinee likelihood function on a uniformly spaced grid of ng points, (1; (q1),91),.--, 
(t (4n, sti) and applies the formulas 


Ao De ah (4) N (q; P21, Z) 


— 4.1 
oe ly (a) N (3 Fits 2) { 


9 


du 


='\2 
Dita (4: -6;) h(a) N (gi; Py2i, 2) 
mn 
Dies (gi) N (gi; Ps. Z) 
This approach is described in Beaton(1987). Gaussian quadrature can be used as an al- 


(4.2) 


‘var (6; | X,Y,I1, £1) = 


ternative which yields similar accuracy with fewer quadrature points, but the Gaussian 
quadrature is still less efficient because it requires that the examinee likelihood function be 
re-evaluated at each E-step, in contrast to the uniformly spaced quadrature where the grid 
of examinee likelihood function evaluations can be performed once before beginning the 
EM iterations. The generalization to a two dimensional grid of nr points for two scales is 
immediate. The numerical quadrature results in Table 1 were computed using a uniformly 
spaced two dimensional grid of 41? points on a square region determined by the values of 
the item specific parameters. Numerical quadrature routines are computationally unfeas:ble 
for applications with more than two scales. 

The Integrated normal approximation referred to in Table 1 is much faster than the other 
computing approximations, but does not produce adequate accuracy for typical applications 
having scales with less than 20 examination items per examinee. The Integrated normal 
approach approximates the contribution to the examinee likelihood function associaced 
with eaci. scale by a normal distribution with mean and variance computed using the grid 
of equally spaced quadrature points, 


bes Yt, gil; (93) 
BOS Gel’ a 
a: itn (Gi — Hy)? (a) 
Sai Lj (91) 


Because the examinee likelihood function factors, the examinee likelihood function is ap- 


7; oe Cees (4.4) 


proximated by a multivariate normal distribution with a diagonal variance-covariance ma- 
trix, 1(@) = N(,D) with w = (11,...,up)’, and D = diag(7,...,7)). The examinee 
posterior mean and variance-covariance matrix in each E-step are computed using Bayesian 


normal theory calculations, Box and Tiao(1973), 
-1 
var (0 | X,Y, I, 5) (71+ D"') =V, (4.5) 


10 
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G=V(FTriz+D"p) . (4.6) 


Mislevy, Johnson, and Muraki(1992), and Mislevy(1990) describe this approach. The In- 
tegrated normal approximation produces more accurate results than the common modal 
approximation in (3.1) and (3.2), and the Integrated normal approximation is much faster 
because the moments for approximating the examinee likelihood function can be evaluated 


once before the EM iterations begin. 
5 Relative speed of the different methods 


Asymptotic corrections to the normal approximation have accuracy close to that of ex- 
tensive product quadrature rules and the corrections are applicable when there are many 
scales. Timing runs show that the evaluation of the partial derivatives of the examinee log 
posterior density for all scales requires about the same amount of time as the evaluation of 
the examinee log posterior density at two distinct values of 8. A Newton-Raphson itera- 
tion thus requires the equivalent of two posterior density evaluations, and likewise for the 
computation of the correction formulas. Three Newton-Raphson iterations yield sufficient 
accuracy for @ resulting in an algorithm requiring the equivalent of eight posterior density 
evaluations. The evaluation of the examinee log posterior density grows linearly in the 
number of examination items (and scales) producing a very efficient algorithm for problems 
involving multi-dimensional 6,. 

The Integrated normal algorithm is approximately ten times as fast as the numerical 
quadrature algorithm i: two dimensions. The Integrated normal algorithm is approximately 
five times as fast as the asymptotic correction algorithm regardless of dimension. Applica- 
tions of the asymptotic correction algorithm typically require one to four hours on a Sun 


SPARC Station IPX. 
6 Overview of the MGROUP program 


MGROUP is a FORTRAN77 program with many subroutines that are called by two primary 
controlling routines called PHASE] and PHASE2. Mislevy and Shechan(1987) provides 


1] 


documentation for earlier versions of the program. The PHASE] routine calls subroutines 
that input data and perform calculations that can be done once before beginning the EM 
iterations. The PHASE2 routine directs the EM iterations, ani then executes subroutines 
to report on the converged values. If requested, PHASE2 also executes subroutines to 
compute standard errors for I and generated imputed values for each examinee. The 
code for PHASE] and PHASE2 are given in Appendices B and C. The most important 
computational steps are briefly outlined in this section; the outlined steps are highlighted 
in the code for easy reference. 

There are currently three versions of the MGROUP program that implement the dif- 
ferent methods of approximating the E-step expectations: BGROUP, which implements 
the numerical quedrature approach in (4.1) and (4.2), NGROUP which implements the 
Integratcd normal approximation in (4.6) and (4.5), and CGROUP which implements the 
asymptotic correction formulas in (3.3) and (3.4). The different versions of the MGROUP 
program share most subroutines. Subroutines that are specialized for different versions of 
MGROUP are indicated by inserting a ’B’, °C’, or 'N’ in the subroutine name, for example, 
BPHASE2.f, CPHASE2.f, NPHASE2.f. The function of the specialized subroutines is the 
same for each version of MGROUP, but the data structures and output are modified to 
accomadate the differing requirements of the E-step calculations. Subroutines displayed in 


the appendices are for CGROUP v.2.2, except where explicitly indicated otherwise. 


PHASE1 


1. Input the item parameters, scale definitions, conditioning variable definitions, and the 


group definitions (if requested). 


2. Input the item response data and the conditioning variables (READDF/READMF). 
For each examinee, evaluate the equally spaced grid of likelihood values described 
in Section 4. Use this grid of points to obtain (4.3) and (4.4). Output a work file 


containing the conditioning variables for use in the E-step. The BGROUP version also 


12 
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outputs the grid of likelihood values for each examinee in the work file. The examinee 
likelihood function is normed to have a maximum value of one. The NGROUP version 
appends only (4.3) and (4.4) to the conditioning variables. The CGROUP version 
appends (4.3) and (4.4) for use in creating starting values to find the mode in (3.2) for 
each examinee in the E-step. The item responses are included in the CGROUP work 
file so that the derivatives of the examinee log likelihood function can be computed. 
The CGROUP version also includes the maximum value of the (un-normed) examinee 
log likelihood function for use in norming the approximate log likelihood function for 


IF and J. 
3. Compute (X’X)7’ for use in the M-step of the EM algorithm. 


4. Input starting values for F and XS or compute them using (4.3) regressed on X using | 


summary statistics accumulated during steps 2 and 3. 


PHASE2 
1. EM iterations. 
(a) Execute the subroutine ESTEP which computes the examinee posterior mean 


and variance and accumulates the required summaries of the moments. 


(b) Execute the subroutine MSTEP which is very simple because (X’X)~' is already 


computed. 


(c) Check the covergence criteria. Two criteria are checked: 1) the largest change 
in I after dividing by the standard deviations of X, 2) the change in the log 


likelihood function. 


2. Produce printed summaries of converged values and write I and Z to a file if re- 


quested. 


3. Produce standard errors for I if requested. 
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4. Generate multiple imputations if requested. The imputations are produced by the 
subroutine WRITSF. This subroutine uses the same code as in the ESTEP to com- 
pute an approximate posterior mean and variance for each examinee; the subroutine 
then generates imputations from a normal random variable with the computed mean 
and variance. Changes to the code of the E-step require corresponding changes in 


WRITSF. 
7 Numerical integration: BGROUP E-step 


The code implementing the steps outlined below is given iz Appendix D. 


1. Read the data from the work file created by READDF/READMF in PHASE1. The 


data are read one examinee at a time and steps 2-8 are completed for each examinee. 


2. Compute the prior mean using I; from the previous iteration and the conditioning 


variables from the PHASE1 work file. 


3. Evaluate the examinee posterior density given in the integrand of (1.4) at each of 
the uniformly spaced quadrature points created in PHASE1. The density is only 


determined up to a norming factor at this point. 


(a) Compute the normal prior density for the first scale. 


(b) If there are two scales, then compute the conditional prior distribution for the 
second scale conditional on the quadrature value of the first scale for each quadra- 


ture point in the second scale. 


(c) When there are two scales, multiply 3a) and 3b) to obtain the prior distribution. 
Multiply the prior by the likelihood function for each scale at the corresponding 
quadrature points. The likelihood function is obtained from the PHASE] work 
file. 


4. Accumulate the first and second moments for (4.1) and (4.2). These moments are 


only determined up to a norming constant. 
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5. Compute the norming factor by summing the posterior density evaluated at each 


quadrature point. 


6. Accumulate the examinee contribution to the log likelihood function in (1.4) approx- 
imating the integral with the uniformly spaced quadrature. Note that I; and Z;, 
which are the estimates from the previous iteration, are used so that the likelihood 


function evaluations always lag the estimates by one iteration. 
7. Center the posterior variance calculations and divide by the norming constant. 


8. Accumulate the posterior means and variances for (1.7) and (1.8). Note that I’, and 


not I’;41 is used in the centering of (1.8). 


8 Integrated normal approximation: NGROUP E-step 


The code implementing the steps outlined below is given in Appendix E. 


1. Read the data from the work file created by READDF/READMF in PHASE]. The 


data are read one examinee at a time and steps 2-6 are completed for each examinee. 


2. Compute the prior mean using I, from the previous iteration and the conditioning 


variables from the PHASE1 work file. 
3. Compute the posterior variance using (4.5). 
4. Compute the posterior mean using (4.6). 


5. Accumulate the posterior means and variances for (1.7) and (1.8). Note that I’, and 


not I”;41 is used in the centering of (1.8). 


6. Accumulate the examinee contribution to the log likelihood function in (1.4). The 
Integrated normal approximation in Section 4 is equivalent to a measurement er- 


ror model with “observation” yz; distributed as N(@,,.D,), and 6, distributed as 


9 


N(I'2;,). Thus, uw; ~ N(I’a;, + D,), and the contribution to the likelihood 
function in (1.4) with (1.1) approximated by this normal density at (I;, 3) is 


F log (27) - 5 log |X: + Di| - 5 (mu; — Ti2i) (21+ Di)! (u;—Piai) . (8.1) 


"Note that (p/2) log (27) is added to (8.1) to obtain consistency with the BGROUP 


likelihood calculation which has the maximum of the examinee likelihood function set 
to one. Also note that I’, and X',, which are the estimates from the previous iteration, 
are used so that the likelihood function evaluations always lag the estimates by one 


iteration. 


Asymptotic corrections: CGROUP E-step 


The code implementing the steps outlined below is given in Appendix F. 


. Read the data from the work file created by READDF/READMF in PHASE1. The 


data are read one examinee at a time and steps 2-8 are completed for each examinee. 


. Compute the prior mean using I, from the previous iteration and the conditioning 


variables from the PHASE1 work file. 


. Complete the NGROUP steps 3 and 4 to produce starting values for finding 0. 


. Make a copy of the NGROUP estimates for use if the asymptotic corrections fail. The 


posterior standard deviations from NGROUP are also used for bounding the amount 
of change allowed in the Newton-Raphson iterations and the asymptotic corrections. 
The current value of the step size for scale J stored in XSTEP(J) is one posterior 
standard deviation. The maximum step size can be changed by modifying the value 


of XNR1 in the subroutine GETGLO. 


. Perform Newton-Raphson iterations to find 8. 


(a) Compute the first and second derivatives of the examinee log likelihood function 


and the log posterior distribution. 
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t ata) 
LO 


(b) Compute the Cholesky decomposition of the matrix of second derivatives of the 


log posterior distribution and use it to compute the Newton-Raphson update. 
The Cholesky subroutine DCHOL returns an appropriately chosen “square root” 
matrix when the matrix of second derivatives is not positive definite; the algo- 
rithm employed is described in Kennedy and Gentleman(1980), p. 442-450, with 


bounds obtained from the NGROUP posterior variance approximation. 


(c) Update the Newton-Raphson approximation to the mode limitirg the change in 


scale J to size XSTEP(J), and check for convergence. The maximum number of 
iterations allowed is MAXIT=6 and the current convergence criteria is XNR2=0.001; 


these values are set in the subroutine GETGLO. 


6. Compute the asymptotic corrections in Section 3. 


(a) Compute G by inverting the second partial derivatives of the log posterior distri- 


bution computed at the last Newton-Raphson iteration. Also compute the third 
and fourth derivatives of the log posterior (which equal to the derivatives of the 


log likelihood) at the converged value of 8. 


(b) Compute the asymptotic corrections to @ in (3.3). 


(c) If the difference between the corrected approximation to the posterior mean and 


the NGROUP approximation to the posterior mean exceeds XNR4*XSTEP, then set 
the corrected value equal to the NGROUP approximation plus the appropriately 
signed value of XNR4*XSTEP. The current value of XNR4 is 1.0 which is set in the 
subroutine GETGLO. 


(d) Determine the bound on the largest change allowed for the posterior variance 


for each scale. Only the diagonal elements of the variance matrix approxima- 
tions are checked. A best. preliminary approximation to each variance is formed 
as follows. Each variance approximation in G is divided by the corresponding 
NGROUP approximation to the variance and if this ratio is greater than XNRS, 


then the NGROUP approximation is used for bounding the asymptotic correc- 


if 


a5 
C5 


tions, and if the ratio is less than XNR5, a weighted average of these two variance 
approximations is used (where the weights used are based on several examples). 
The current value of XNRS5 is 1.5 which is set in subroutine GETGLO. 


(e) Compute the asymptotic corrections to G given in (3.4). 


(f) Bound the cize of the asymptotic corrections to the variance. If the ratio of the 
corrected variance to the best preliminary approximation exceeds XNR5, the cor- 
rected variance is set to this upper boundary value. If the ratio of the corrected 
variance to the best preliminary approximation is less than XNR6, then the cor- 
rected variance is set to this lower boundary. The current value of XNR6 is 0.25 
which is set in subroutine GETGLO. 


(g) If any asymptotic correction to the posterior variances fails, then use the corre- 
lations from the NGROUP approximation along with the (bounded) corrected 
variances to form a combined positive definite approximation to the posterior 


variance. 


7, Accumulate the posterior means and variances for (1.7) and (1.8). Note that I, and 


not I,4; is used in the centering of (1.8). 
8. Compute the contribution of each examinee to the log likelihood function. 
(a) If 6 could not be found using Newton-Raphson iterations, then apply the calcu- 
lations from NGROUP step 6 for the log likelihood function. 


(b) Calculate the approximate likelihood function using the usual normal modal 


approximation that omits all terms of order N~}/? anc smaller, 


log {1 (6) } ~ max log {1 (6)}] + 5 (log IG] - log |) - 
1 


; (6- Iiz,) 5,1 (8- Pia.) . (9.1) 


The first term in (9.1) is the log of (1.1) evaluated at the mode of the posterior 


distribution. It is obtained from the final Newton-Raphson itcration. ‘The second 
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term in (9.1) is the maximum of the log of the likelihood in (1.1). This value 
is computed in READDF/READMF and placed in the PHASE] work file. It 
differs from the first term because @ maximizes the log posterior distribution 
and not the log likelihood function. Note that I; and not I;4; is used in the 
centering of (1.8). 


A Appendix 


This appendix sketches the derivation of the second order correction for the variance in 
(3.5) when there is only one scale. Generalizations of the derivation to the multi-scale case 
are discussed briefly. The extensive calculations required for the expansions were completed 
using the symbolic algebra program Mathematica. 

Following the derivation of Mosteller and Wallace(1964), p. 153, make the change of 
variable u = K}/? (41 - 61) , and modify the definition of h by the factor K~! for convenience 
in tracking the order of the terms, h = —K™!log[f (y | @) N(@;2'I';, Z;)]. The Taylor 


expansion of h about 6; is 


Gyw . A® us . AM“ us aie 
2 6K1/2 24K , 


Kh(6:) = Kh (6) + 


recalling that Gj, is the inverse of the second derivative of h evaluated at 6. The second 


non-central posterior moment can thus be written as 


~Kh(6;) 9249 
2\ _ fe 1271 
E (0) ~ fenKhG)d0, 


A y3 AM us 


(27Gy,)7'/? f exp (-}G;)u?) exp (-Sos -tr-" : (6 +2K-/?6,u + K-v?) du 
(27G}}) a fexp (-3G;}u2) exp (-s - aie ~ -:) du 
(A.1) 
The integrals in the numerator and denominator of (A.1) can be viewed as the expectations 
of complicated functions of a normal random variable U with mean 0 and variance G)). 


Expanding 


4 
exp Graz aie ae ) (A.2) 


as a polynomial in u and computing the resulting higher moments of U term wise produces 
series representations for the numerator and denominator in (A.1). The Taylor series ex- 
pansion of a ratio is then applied to the series representations for (A.1) and terms smaller 


than K~? are deleted yielding 
e (9) = OF + KT! (Gi - O,Gi, hy") +k? (3) Gi; [1561 (i”)" 
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iz . 3 a a « Py ror 
-156,63, (A)? ~ 6h” + 16 GnhOAM — 36, i (A.3) 


Similar calculations for E (@;) produce the second order approximation 
es = 1 a 3 | . 3 
E((i) = 6—-K 1 (5) G2h®) +K~? (3) G3, [- 15G3, (Ai) 


Applying the identity var (@,) = E (6?) —(E (6:))? to (A.3) and (A.4), deleting terms smaller 
than K a8, and redefining h without K~! produces (3.5). (Note: The refined correction in 


(A.4) has not been implemented because (3.3) has performed well in practice and (A.4) 
requires triple summations in the correction formulas in the multi-scale case and the com- 
putation of fifth derivatives). 

The multi-scale case proceeds similarly with transformed variables U;,...,Up that have 
a multivariate normal distribution with mean zero and variance-covariance matrix with 
elements G;,. There is an exponential term of the higher powers of u like (A.2) for each 
uy in the multi-scale case, but this is still a tremendous simplification of the general second 
order expansion because there are no mixed (u;, u,) terms of order three or higher in the 
multivariate Taylor expansion of h. The remaining complication in the multi-scale case 
is that the calculation of the higher univariate normal distribution moments required for 
the single scale case are replaced by numerous high order multivariate normal distribution 


moment calculations in the multi-scale case. 
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B Appendix Subroutine PHASE1 


Ceeee Sse ee SSS SESS SSSSSSSSSSSSSSSSSSSSSSESSSSSSSSSSESSSSSSSERSSESSSS 


C SUBROUTINE PHASE1 


C PURPOSE: SET DEFAULTS AND/OR READ INPUT FILES FOR STARTING 

C ESTIMATION OF GAMMA 

Cc 

CSS Se eee ee SS SSS BSSSSSSSSSSSSESSSSESSSESSSSSSSSRSSSSSSSSSSSSSSSESSSSSSE 
SUBROUTINE PHASE1(#) 
INCLUDE ’BIGMAT.FOR’ 
INCLUDE ’GLOBAL.FOR’ 

Cc 


C PHASE] 1. 


Cc 

C##NT 

C CHANGE TO NUMERIC CALL 

C 
INTEGER*4 JTIME(3) 

Cc 

Cc 
ICYCL = 0 
CALL ALLOC(KINAM,8,NITM) 
CALL ALLOC(KISLO,8,NITM) 
CALL ALLOC(KITHR,8,NITM) 
CALL ALLOC(KIGUE,8,NITM) 
CALL ALLOC(KIUSE,4,NITM) 
CALL ALLOC(KILOC,4,NITM) 
CALL ALLOC(KILNK ,4,NITM) 
CALL ALLOC(KIKEY ,4,NITM) 
CALL ALLOC(KIOMT,4,NITM) 
CALL ALLOC(KINPR,4,NITM) 


SPECIAL TO CGROUP VERSION OF THE PROGRAM 


aaa 


CALL ALLOC(KOSUB,4,NITM) 


#NT ALLOCATE POLYTOMOUS ITEM PARAMETERS 
CHANGED TO MEET AL’S NEW CALLING SEQUENCE 


aqaaqaa 


CALL ALLOC(KICAT,8,NITM*MAXCAT) 
CALL ALLOC(KIALT,4,NITM) 

CALL ALLOC(KINCA,4,NITM) 

CALL ALLOC(KTNAM,8,NTEST) 

CALL ALLOC(KTLEN ,4,NTEST) 


CALL READIF(CM(KINAM) ,BM(KISLO) ,BM(KITHR) ,BM(KIGUE) , 
$ BM(KICAT) , IM(KIUSE) , IMCKILOC) , IM(KILNK) , 
$ IM(KIALT) , IM(KINCA) , IMC(KIKEY) , IMCKIOMT) , 
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Q 


Cc 


$ IM(KINPR) ,CM(KSNAM) , IM(KSLEN) ,SM(KOSUB) , 
$ CM(KTNAM) , IM(KTLEN) , #500) 


IF (INSCAL .GT. 0) THEN 
IBTEMP = IBASE 
CALL ALLOC(IWRK1 ,4,NTOTI) 
CALL GETSCA(CM\KSNAM) , IM(KILOC) , IM(KILNK) , IM(IWRK1) , #500) 
IBASE = IBTEMP 
END IF 


#NT CHANGED TO NEW PRINT ROUTINE PER A. ROGERS 


NLEFT = 0 

CALL PRTSCA(CM(KSNAM) , IM(KSLEN) , 

$ CM(KINAM) ,BM(KISLO) ,BM(KITHR) , BM(KICUE) ,BM(KICAT) , 
$ IM(KIUSE) , IM(KILNK) , IM(KILOC) , IM(KIALT) , IM(KINCA) , 
$ IM(KIKEY) , IM(KIOMT) , IM(KINPR) ) 


MGRP = MAX(NGRP, 1) 
MCYCLE = MAX(1,NCYCLE+1) 
CALL ALLOC(KCNAM,8,NRANK) 
CALL ALLOC(KCUSE,4,NRANK) 
CALL ALLOC(KCLOC,4,NRANK) 
CALL ALLOC(KCLAB,8,5*NRANK) 
CALL ALLOC(KGNAM,8,MGRP) 
CALL ALLOC(KLNAM,8 ,MAXLEV*MGRP) 
CALL ALLOC(KGLEV 4 ,MGRP) 
IF (CFDISP .EQ. 1 .OR. CFDISP .EQ. 3) THEN 
CALL READCF(CM(KCNAM) , IM(KCUSE) , IM(KCLOC) ,CM(KCLAB) , 
$ CM(KGNAM) , CM(KLNAM) , IM(KGLEV) , #500) 
ELSE 
IF (INCOND .GT. 0) THEN 
CALL GETCON(CM(KCNAM) , IM(KCUSE) , IM(KCLOC) ,CM(KCLAB) , #500) 
ELSE 
CALL SETCON(CM(KCNAM) , IM(KCUSE) , IM(KCLOC) ,CM(KCLAB) ) 
END IF 
IF (NGRP .GT. 0 .AND. INGRP .GT. 0) THEN 
CALL GETGRP(CM(KGNAM) ,CM(KLNAM) , IM(KGLEV) , *500) 
ELSE 
CALL SETGRP(CM(KGNAM) ,CM(KLNAM) , IM(KGLEV) ) 
END IF 
END IF 


C PHASE1 1. END 


Cc 


CALL ALLOC(KQPT,8,NQPT) 

CALL ALLOC(KTDT,8 , NRANK*NRANK) 
CALL ALLOC(KCOV,8,NSCALE*NSCALE) 
CALL ALLOC(KTDM,8 ,NRANK*NSCALE) 
CALL ALLOC(KPMN,8 ,NSCALE) 
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aagaaaaa 


aa 


Qa CP ae 


Q 


CALL ALLOC(KPVAR ,8,NSCALE) 

CALL ALLOC(KQMN ,8,NSCALE) 

CALL ALLOC(KQVAR ,8,NSCALE) 

CALL ALLOC(KTMU,8,NRANK) 

CALL ALLOC(KTSD, 8, NRANK) 

CALL ALLOC(KLK 8 ,NQPT*NSCALE) 

CALL ALLOC(KGUN ,8, TOTLEV) 

CALL ALLOC(KGWN,8, TOTLEV) 

CALL ALLOC(KGMU, 8, NSCALE*TOTLEV*MCYCLE) 
CALL ALLOC(KGSD,8,NSCALE*TOTLEV*MCYCLE) 
CALL ALLOC(KTVEC ,4,NRANK) 

CALL ALLOC(KGVEC ,2,NTOTG) 

CALL ALLOC(KILO,4,NSCALE) 

CALL ALLOC(KIHI ,4,NSCALE) 


ADDED BY NT 9/29/91 


VARIABLES FOR THE STDERR ROUTINE THAT ARE ACCUMULATED 
IN READDF AND READMF 


CALL ALLOC(KTAU,8,NSCALE) 
CALL ALLOC(KTT,8,NRANK*NRANK) 


VARIABLE ADDED TO ALLOW DIRECT COMPUTATION OF 
DERIVATIVES OF THE LIKELIHOOD FUNCTION. 

CALL ALLOC(KXITEM,4,NITM) 

IBTEMP = IBASE 

CALL ALLOC(KCOR, 8, NRANK*NRANK) 

## CHANGED 5/4/92 TO CONFORM TO NEW POLYTOMOUS CATEGORY USAGE 
CALL ALLOC(KLOGP ,8 , NQPT*NITM* (MAXCAT+2) ) 

CALL ALLOC(KLKI ,8,NQPT) 

CALL ALLOC(KQLAB ,8,NQPT) 

CALL ALLOC(KXVEC,2,NITM) 

CALL ALLOC(KNTRY ,4,NSCALE) 


## ADDED 10/1/92 FOR ITEM RESPONSE FREQUENCIES 


PHASE] 2. 


CALL ALLOC(KFREQ,8,'%ITM* (MAXCAT+3) ) 
CALL ALLOC(KPCT,8,NITM*(MAXCAT+1) ) 


CALL SETQPT(BM(KQPT ) ,BM(KISLO) ,BM(KITHR) ,BM(KIGUE) ,BM(KICAT), 
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$ IM(KIUSE) , IM(KIALT) , IM(KINCA) ,BM(KLOGP) ,BM(KLKI )) 
IF (NDIAG .GE. 3) THEN 
CALL PRTMAT(BM(KLOGP) ,NQPT,CM(KQLAB) ,NQPT,CM(KINAM) ,NITM, 
$ *TABLE OF LOG(CORRECT) PROBABILITIES BY ITEM’) 
CALL PRTMAT(BM(KLOGQ) ,NQPT ,CM(KQLAB) ,NQPT ,CM(KINAM) ,NITM, 
*TABLE OF LOGCINCORRECT) PROBABILITIES BY ITEM’) 
CALL PRTMAT(BM(KLKI ) ,NQPT,CM(KQLAB) ,NQPT,CM(KSNAM) ,1, 


qgqaagaaaana 
“A 


$ * INITIAL LIKELIHOOD ESTIMATES’ ) 
END IF 
IUNIT = 1 
Cc 
IF (DFDISP .EQ. 1) THEN 

CALL GETINP 

CALL ALLOC(KXC,4,NTOTC) 

CALL ALLOC(KXI,4,NTOTI) 

CALL READDF(BM(KQPT ),BM(KTDT ),BM(KCOV ),BM(KTDM ),BM(KPMN ), 
$ BM(KPVAR) ,BM(KQMN ),BM(KQVAR) ,SM(KTVEC) ,BM(KLK ), 
$ IM(KILO ),IM(KIHI ),BM(KLOGP) ,BM(KLKI ), 
$ HM(KXVEC) ,SM(KXC  ),IM(KXI ),IM(KNTRY), 
$ IM(KCUSE) , IM(KCLOC) , IM(KIUSE) , IM(KILOC) , IM(KILNK) , 
$ IM(KINCA) , IM(KIKEY) , IMC(KIOMT) , IM(KINPR) ,CM(KSNAM) , 
$ HM(KGVEC) , IM(KGLEV) ,BM(KGUN ) .BM(KGWN ),BM(KGMU ), 
$ RM(KGSD ) ,BMC(KFREQ) ,BM(KTAU ),BM(KTT ),SM(KXITEM), 
$ SM(KOSUB) , #500) 

ELSE IF (MFDISP .EQ. 1) THEN 

CALL READMF(BM(KQPT ) ,BM(KTDT ),BM(KCOV ),BM(KTDM ),BM(KPMN ), 
$ BM(KPVAR) ,BM(KQMN ) ,BM(KQVAR) ,SM(KTVEC) ,BM(KLK ), 
$ IM(KILO ),IM(KIHI ),BM(KLOGP) ,BM(KLKI ), 
$ HM(KXVEC) , IMC(KNTRY) , IM(KCUSE) , IM(KILNK) , IM(KINCA) , 
$ CM(KSNAM) ,HM(KGVEC) , IM(KGLEV) ,BM(KGUN ),BM(KGWN ), 
$ BM(KGMU ) ,BM(KGSD ) ,BM(KFREQ), 
$ BM(KTAU ),BM(KTT ),SM(KXITEM) ,SM(KOSUB) , #500) 

END IF 
Cc 
C PHASE1 2. END 
Cc 
Cc 
c PRINT ITEM FREQUENCY STATISTICS ADDED 10/1/92 
Cc 

CALL STATS(CM(KTNAM) ,CM(KINAM) , IM(KILNK) , IM(KIUSE) , IM(KILOC) , 

$ IM(KINCA) ,BM(KFREQ) , BM(KPCT) , #500) 
Cc 
Cc 
C PHASE1_ 3. 
Cc 

IBASE = IBTEMP 

Cc 


qd 
ae 


KTT IS THE UNWEIGHTED DESIGN MATRIX USED FOR STANDARD ERROR 
AND IMPUTATION OF GAMMA 


aaa 


IF (NRANK .GT. 1) THEN 
CALL PRIDES(BM(KTDT ) ,CM(KCNAM) ,BM(KCOR ) ,BM(KTMU ),BM(KTSD )) 
CALL INVDES(BM(KTDT ) ,BM(KCOR ) , IM(KCUSE) , #500) 
CALL PRTCON(CM(KCNAM) , IMCKCUSE) , IM(KCLOC) ,BM(KTMU ) ,BM(KTSD ), 
$ CM(KCLAB) ) 


ADDED BY NT 9/29/91 


aaa 


CALL INVDES(BM(KTT ),BM(KCOR ) ,IM(KCUSE) , *500) 


ELSE 
BM(KTDT) = 1D0/BM(KTDT) 


ADDED BY NT 9/29/91 


aaa 


BM(KTT) =1D0/BM(KTT) 
END IF 


c 
C PHASE1 3. END 
¢ 
C 


IF (NGRP .GT. 1) THEN 
CALL PRTGRP(CM(KGNAM) ,CM(KLNAM) , IM(KGLEV) ,CM(KSNAM) , 
$ BM(KGUN ),BM(KGWN ),BM(KGMU ),BM(KGSD )) 
END IF 


IF (CFDISP .EQ. 2 .OR. CFDISP .EQ. 3) 


$ CALL WRITCF(CM(KCNAM) , IM(KCUSE) , IM(KCLOC) ,CM(KCLAB) , 
$ CM(KGNAM) , CM(KLNAM) , IM(KGLEV) , #500) 


PHASE1 4. 


a Q ee 


CALL ALLOC(KGAM,8,NRANK*NSCALE) 
CALL ALLOC(KSIG,8,NSCALE*NSCALE) 
CALL ALLOC(KFIXV,8,NSCALE) 
IF (GFDISP .EQ. 1 .OR. GFDISP .EQ. 3) THEN 
CALL READGF(BM( KGAM) ,BM( KSIG) ,BM(KFIXV) ,*500) 
ELSE 
CALL SETGAM(BM( KGAM) ,BM( KSIG),BM( KTDT) ,BM( KTDM), 
$ BM( KCOV) ,BM(KQVAR) , BM(KFIXV) ) 
END IF 


aQ 


C PHASE1 4. END 
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NLEFT = 0 
CALL PRTGAM(BM( KGAM) ,BM( KSIG) ,BM( KCOV) ,BM(KPVAR), 
$ CM(KCNAM) ,CM(KSNAM) ) 


C##NT 
Cc CALL GETTIMCIHR, IMIN, ISEC, IHUN) 
CALL ITIME( JTIME) 
IBHIGH = 4*(IBMAX~-1) 
KALLOC = 4*#IALLOC 
WRITE(6,9200) JTIME(1) , JTIME(2) , JTIME(3) , IBHIGH, KALLOC 
Cc WRITE(6,9200) TCHAR, IBHIGH ,KALLOC 


IBMAX = IBASE 
RETURN 
Cc 
500 WRITE(6,9000) 
RETURN 1 
9000 FORMAT(’ *** EXECUTION TERMINATED IN SUBROUTINE PHASE1 ««’) 
9200 FORMAT(’OPHASE 1 COMPLETED AT ’,12.2,’:’,12.2,’:’,12.2, 


$ » USING ’,I8,’ OF ’,I8,’ BYTES OF WORKSPACE’ ) 
C9200 FORMAT(’OPHASE 1 COMPLETED AT ’,A8, 
Cc $ ’ USING ’,18,’ OF ’,I8,’ BYTES OF WORKSPACE’ ) 
C#NT THRU HERE 
END 
al 


ao 
Co 


Q 


aqgaagqagqaaana 


Appendix Subroutine PHASE2 


VERSION OF PHASE2.F CREATED FOR CORRECTED ESTEP 


1/9/92 NEAL THOMAS 


SUBROUTINE PHASE2(*) 
IMPLICIT REAL*8 (A-H,0-Z) 
INCLUDE ’BIGMAT.FOR’ 
INCLUDE ’GLOBAL.FOR’ 


C##NT 


C 
C 


Qaaaaaa 
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CHANGE TO NUMERIC CALL 


INTEGER*4 IHR, IMIN, ISEC, IHUN 
INTEGER*4 JTIME(3) 


CHARACTER*30 HEDGAM, HEDSIG 
CHARACTER*80 HEDCHG 


##NT VARIABLES FOR REPORTING ON THE PERFORMANCE OF THE 
CORRECTION TO THE NORMAL APPROXIMATION TO THE ESTEP. 
PASSED TO PHASE2 SO THAT THEY CAN BE REPORTED AT THE 
‘END OF THE PRINTING 


INTEGER ISUBJ, IFAIL, INOTPD, IFCORM, IFCORV 
REAL*8 AVEIT 


DATA HEDGAM /’ ACCELERATORS FOR GAMMA tf 

DATA HEDSIG /’ ACCELERATORS FOR SIGMA +} 

DATA HEDCHG 

$ /’CHANGES IN STANDARDIZED GAMMA FROM PREVIOUS CYCLES FOR SCALE:’/ 


SDIM = DFLOAT(NSCALE) 

GDIM = ODO 

DO 10 I=1,NRANK 

IF (IM(KCUSE+I-1) .EQ. 1) GDIM = GDIM+SDIM 

CONTINUE 

IBTEMP = IBASE 

MCYCL = MAX(1,NCYCLE+1) 

CALL ALLOC(MGAM,8, NRANK*NSCALE*MAX (MCYCL ,NSIMP) ) 
CALL ALLOC(MSIG,8,NSCALE*NSCALE*MCYCL) 

CALL ALLOC(KGACC,8 , NRANK*NSCALE) 

CALL ALLOC(KSACC,8 , NSCALE) 

CALL ALLOC(KCYCL,8 ,MCYCL) 

CALL SETACC(BM(KGAM) ,BM(MGAM) ,BM(KGACC) , 
$ BM(KSIG) ,BM(MSIG) ,BM(KSACC) , 
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$ CM(KCYCL) ) 


1 LKTOT FOR EACH CYLCLE, CHANGED BY NEAL THOMAS 


Qagaa 


CALL ALLOC(KLKT,8,MCYCL) 

CALL ALLOC(KSVAR,8,NSCALE) 

CALL ALLOC(KCOVI ,8 ,NSCALE*NSCALE) 
CALL ALLOC(KCOR,8,NSCALE*NSCALE) 
CALL ALLOC(KINDEX ,4,NSCALE) 

CALL ALLOC(KID,4,NSCALE) 


CALL ALLOC(KD,8,NSCALE) 

CALL ALLOC(KGAMCH,8, NRANK*NSCALE*NLAG) 
CALL ALLOC(KBIGG,8,NLAG) 

CALL ALLOC(KSUMG ,8,NLAG) 


##NT 
ALLOCATE SPACE FOR DERIVATIVES OF THE LOGLIKELIHOOD AND POSTERIOR 


aQqaaaa 


CALL ALLCC(KDERIV,8,4*NSCALE) 
CALL ALLOC(KGSP ,8,NSCALE*NSCALE) 
CALL ALLOC(KBACK ,8,NSCALE) 

CALL ALLOC(KCB,8,NSCALE*NSCALE) 
CALL ALLNC(KMODE,8,NSCALE) 


Q 


##ADDITIONAL WORK SPACE 


CALL ALLOC(KXMI,8,NSCALE) 

CALL ALLOC(KXSTEP, 8 .NSCALE) 

CALL ALLOC(KSINV,8,NSCALE*NSCALE) 
CALL ALLOC(KWMEAN ,8,NSCALE) 


C 
c ##NT = WORK SPACE FOR LABEL FOR PRTSUM 


Cc 
CALL ALLOC(KRLAB ,8,MCYCL) 
C 
C 
CesesseesesrserestrretessseassssssresssesssEsarsnsssssssseeseecss= 


C NOW BEGIN THE E-M CYCLE 

Ceesserecsssreeereeeessrrsrrstsssrassiarsssesesesssssesssesssssse= 
IRS = NRANK*NSCALE 

ISS = NSCALE*NSCALE 

JGAM = MGAM 

JGAM1 = JGAM-IRS 

JGAM2 = JGAMi-IRS 

JSIG = MSIG 

JSIGi1 = JSIG-ISS 

JSIG2 = JSIG1-ISS 


C CHANGED BY NEAL THOMAS SO THAT ALWAYS ONE LKTOT PER ITERATION 
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, 


qe 


JLKT = KLKT+1 
JGMU = KGMU+NSCALE*TOTLEV 
JCYCL = KCYCL+MCYCL-1 


Cc 
C ##NT ASSIGN FLAG INDICATING THE RUN IS A RESTART 
Cc 
IRST=0 
IF( (GFDISP.EQ.1).OR.(GFDISP.EQ.3) )IRST=1 
C 
c ##NT 
c LIKELIHOOD CALCULATIONS LAG PARAMETER ESTIMATES BY ONE ITERATION 
Cc SO THAT THE FIRST ENTRY IN THE LIKELIHOOD SUMMARY IS NOT 
Cc ASSIGNED UNLESS THE RUN IS A RESTART 
Cc 
IF (IRST.EQ.1) THEN 
BM(KLKT)=XLF 
XLNORM=XLI 
ELSE 
XLNORM= 
END IF 
C 
c 
IF (NCYCLE .EQ. 0) GO TO 300 
Cc 
Cc 
C PHASE? 1. 
C 
100 ICYCL = ICYCL+1 
JCYCL = JCYCL-1 
Cc 


C PHASE2 1(A) 


CALL ESTEP(BM(KCOV ),BM(KTDM ),BM(KPMN ), 
BM(KXMI ),BM(KQMN ),BM(KQVAR ),SM(KTVEC ), 
BM(JLKT ), 
BM(JGAM ),BM(JSIG ),BM(KXSTEP), 
BM(KCOVI ),BM(KSINV ), 
BM (KWMEAN) , BM(KD ) ,HM(KGVEC ), 
BM(JGMU ),BM(KGWN ),BM(KISLO ),BM(KITHR ), 
BM(KIGUE ),BM(KICAT ),IM(KILNK ),IM(KINCA ), 
SM(KXITEM) ,BM(KDERIV) , 
BM(KGSP ),BM(KBACK ),BM(KCB ),BM(KMODE ), 
ISUBJ, IFAIL, INOTPD, IFCORM, IFCORV, AVEIT) 


PRAPARAAAHAHAA SH 


PHASE2 1(A) END 


30 


UPDATE POINTERS 
JGAM2 = JGAM1 

JGAM1 = JGAM 

JGAM = JGAM+IRS 

JSIG2 = JSIG1 

JSIG1 = JSIG 

JSIG = JSIG+ISS 


C PHASE2 1(B) 


Cc 

CALL MSTEP(BM(KTDT ),BM(KCOV ),BM(KTDM ),BM(KFIXV ), 

$ BM(JGAM ),BM(JSIG ),BM(KSVAR )) 
ACCELERATE GAMMA 

CALL ACCGAM(BM(JGAM ),BM(JGAM1 ),BM(JGAM2 ),BM(KGACC ), 

$ BM(KGAMCH) ,BM(KTSD ),BM(KBIGG ),BM(KSUMG )) 
ACCELERATE SIGMA 

CALL ACCSIG(BM(JSIG ),BM(JSIG1 ),BM(JSIG2 ),BM(KSACC ), 

$ BM(KSVAR ) ,BIGS,SUMS) 


PHASE2 1(B) END 


NLEFT = 0 
CALL PRTGAM(BM(JGAM ),BM(JSIG ),BM(KCOR ),BM(KSVAR ), 
$ CM(KCNAM ) ,CM(KSNAM )) 


IF (NDIAG .GE. 1 .AND. NACC .GT. 0) THEN 
IF(MODC(ICYCL,NACC) .EQ. 0)THEN 
CALL PRTMAT(BM(KGACC) , NRANK ,CM(KCNAM) , NRANK, 


$ CM(KSNAM) , NSCALE , HEDGAM) 
IF CICYCL .EQ. NACC) 
$ CALL PRTMAT(BM(KSACC) ,1, DIAGONAL’ ,1,CM(KSNAM) ,NSCALE, 
$ HEDSIG) 
END IF 
END ‘IF 


SUMS = SUMS/SDIM 
KLAG = MINCICYCL,NLAG) 
IF (NDIAG .GE. 1) THEN 
JGAMCH = KGAMCH 
DO 180 I=1,NSCALE 
WRITE (HEDCHG (63:70) ,8000) CM(KSNAM+I-1) 
CALL PRTMAT(BM( JGAMCH) , NRANK*NSCALE ,CM(KCNAM) , NRANK, 
$ CM(JCYCL) , KLAG ,HEDCHG) 
JGAMCH = JGAMCH+NRANK 
180 CONTINUE 
END IF 
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(Se) 
a 


aaa 


aaa 


qgaqaaqaaaqaagaaaaaa 7 mee 


aagaa 


Cc 
Cc 


Cc 


Cc 


NORM -2*LOGLIKELIHOOD BY ITS INITIAL VALUE 
BM(JLKT)=BM( JLKT) -XLNORM 

##NT CHANGED TO ALSO OUTPUT -2*LOGLIKELIHOOD AND DIAGNOSTICS 
IF (NLEFT .LT. 16) CALL HEDOUT 

WRITE(6,9010) ICYCL, (CM(JCYCL+ILAG~-1) , ILAG=1,KLAG) 
WRITE(6,9020) (BM(KSUMG+ILAG-1) , ILAG=1,KLAG) 

WRITE(5,9030) (BM(KBIGG+ILAG-1) , ILAG=1,KLAG) 


WRITE(6 ,9040) BM(JLKT) 
WRITE(6 ,9050) ISUBJ, IFAIL, INOTPD , IFCORM, IFCORV , AVEIT 


CHANGED BY NEAL THOMAS SO THAT ALWAYS ONE LKTOT PER ITERATION 
JLKT = JLKT+1 


JGMU = JGMU+NSCALE*TOTLEV 
BIGG = DABS(BM(KBIGG) ) 


PHASE2 1(C) 


##NT CHECK CONVERGENCE 
ADDED CRITERIA BASED ON -2*L0G(LIK) 
##NT MODIFIED CONVERGENCE CRITERIA AGAIN 10/22/92 


FIRST ITERATION CONVERGENCE CHECKS HANDLED DIFFERENTLY 
DEPENDING ON WHETHER THE RUN IS A NEW JOB OR RESTART 


NEW RUN (DON’T CHECK LIKELIHOOD CRITERIA) 
NOTE THAT XLNORM IS O UNTIL SET HERE FOR NEW RUNS 


IF( (IRST.EQ.0). AND. CICYCL.EQ.1) )THEN 
XLNORM=BM( JLKT-1) 
XLI=XLNORM 
BM( JLKT-1)=BM(JLKT-1) -XLNORM 
RESTART OR ICYCL GREATER THAN 1 
ELSE IF( (DABS(BM(JLKT-2)-BM(JLKT-1)).GT.CRITL) .AND. 
$ (ICYCL .LT. NCYCLE) ) THEN 


GO TO 100 
END IF 


IF (BIGG .GT. CRIT .AND. ICYCL .LT. NCYCLE) GO TO 100 
PHASE2 1. END 
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aa 


200 CONTINUE 


PHASE2 2. 


##NT UPDATE FINAL LIKELIHOOD FOR EXTERNAL STORAGE 


aaa Q ¢a 


XLF=BM(JLKT-1) 


Qa 


NCYCLE = ICYCL+1 
IF (GFDISP .EQ. 2 .OR. GFDISP .EQ. 3) THEN 
CALL WRITGF(BM(JGAM ),BM(JSIG )) 

END IF 

CALL ALLOC(KFLAG,4,NCYCLE*NRANK) 

CALL PRTSUM(BM(MGAM ),BM(MSIG ),BMC(KLKT ),BM(XTSD ), 
CM(KSNAM ),CM(KCNAM ),SM(KFLAG ), IM(KGLEV ), 
CM(KGNAM ),CM(KLNAM ),BM(KGMU ),CM(KRLAB )) 


AA 


PHASE2 2. END 


#NT 
FINAL VALUES FOR GAMMA AND SIGMA COPIED INTO WORK 
SPACE INITIALLY SET ASIDE FOR THESE VALUES. THIS ALLOWS 
THE ARRAYS HOLDING THE HISTORY OF THESE VARIABLES TO BE 
REUSE? AS WORK SPACE BY WRITSF 


agaadaagaagaaana Q 2? 


DO 210 KK=0,NSCALE*NRANK-1 
BM(KGAM+KK) = BM(JGAM+KK) 
210 CONTINUE 
DO 220 KK=0,NSCALE*NSCALE-1 
BM(KSIG+KK) = BM(JSIG+KK) 
220 CONTINUE 


c 
c 
C PHASE2_ 3. 
c 


IF (SE .EQ. 1) THEN 


Cc 
Cc # NT ON 9/29/91 
Cc 

CALL ALLOC(KGSE,8,NRANK*NSCALE) 
Cc 

CALL STDERR(BM( KTT ),BM(KGAM ),BM( KSIG),BM( KTAU), 

$ BM( KGSE),CM( KCNAM) ,CM( KSNAM) ) 
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C 


END IF 


C PHASE2 3. END 


Cc 
300 
Cc 


CONTINUE 


C PHASE2 4. 


Cc 


aa 


qgqaaqaaaaaa 


PARAARAPAAHAHAA SA GH 


Qe 


#ENT 


QqaQqaa 


IF (SFDISP .EQ. 2) THEN 
CALL ALLOC(KSIGE, 8 ,NSCALE*NSCALE) 
CALL ALLOC(KSIGI,8,NSCALE*NSCALE) 
CALL ALLOC(KBETA,8,NSCALE*NRANK) 
CALL ALLOC(KBIMP ,8,NSCALE*NRANK) 
CALL ALLOC(KTHETA ,8,NSIMP*NSCALE) 


CHANGED BY NT FOR USE IN MODIFIED WRITESF 
MAXRS=MAX (NSCALE , NRANK) 
CALL ALLOC(KZ,8,MAXRS) 


THE VARIABLE TDT IS REUSED AS WORKSPACE FOR 

WRITSF IN ORDER TO AVOID ALLOCATING ANOTHER 

NRANK BY NRANK MATRIX. THE VARIABLE HOLDING 

THE GAMMA’S ESTIMATED FOR EACH CYCLE, MGAM, IS ALSO 
REUSED BY WRITSF. 


CALL WRITSF(BM(KQVAR ),BM(KXMI ),BM(KTAU ), 
SMC(KTVEC ),BM(MGAM ),BM(KSIG ), 
BM(KPMN ),BM(KCOVI ),BM(KCOR ),IM(KINDEX), 
IM(KID ), BM(KZ ) ,BM(KTHETA) ,BM(KGAM ), 
BM(KSIGE ),BM(KSIGI ),BM(KBETA ),BM(KTT ), 
BM(KTDT ),BM(KBIMP ),BM(KSINV ),BM(KQMN ), 
BM(KWMEAN) , BM(KXSTEP) ,BM(KD ),BM(KISLO ), 
BM(KITHR ),BM(KIGUE ) ,BM(KICAT ) ,IM(KILNK ), 
IM(KINCA ) ,SM(KXITEM) ,BM(KDERIV) ,BM(KGSP ), 
BM(KBACK ),BM(KCB ),BM(KMODE ),HM(KGVEC ), 
MAXRS) 

END IF 


PHASE2 4. END 


CALL GETTIM(IHR, IMIN, ISEC, IHUN) 
CALL ITIME(JTIME) 

IBHIGH = 4*(IBMAX-1) 

KALLOC = 4*IALLOC 
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40 


WRITE(6,9200) JTIME(1) , JTIME(2) , JTIME(3) , IBHIGH , KALLOC 
Cc WRITE(6,9200) TCHAR, IBHIGH, KALLOC 

IBASE = IBTEMP 

RETURN 


8000 FORMAT(A8) 
9010 FORMAT(’-SUMMARY OF STANDARDIZED GAMMA CHANGES AFTER CYCLE’ ,14/ 


$ 0 FROM :’,T13,10(4X,A8)) 
9020 FORMAT(’ AVERAGE :’,T13,10F12.5) 
9030 FORMAT(’ LARGEST :’,T13,10F12.5) 
9200 FORMAT(’OPHASE 2 COMPLETED AT ’,12.2,’:’,12.2,’:’,12.2, 


$ » USING ’,18,’ OF ’,18,’ BYTES OF WORKSPACE’ ) 
C9200 FORMAT(’OPHASE 2 COMPLETED AT ’,A8, 
Cc $ » USING ’,18,’ OF ’,1I8,’ BYTES OF WORKSPACE’) 


9040 FORMAT(’0 -2*LOG(LIKELIHOOD)=’ ,F14.5) 

9050 FORMAT(’ONUMBER OF SUBJECTS ’,T60,16,/, 
$ ? SUBJECTS WHERE MODE NOT FOUND’ ,T60,16,/, 
$ ? SUBJECTS WHERE INFORMATION NOT POSITIVE DEFINITE’, 
$T60,16,/,’ SUBJECTS WHERE CORRECTION TO MEAN EXCEEDED BOUND’, 
$T60,16,/,’ SUBJECTS WHERE CORRECTION TO VARIANCE’, 
$ » EXCEEDED BOUND’ ,T60,16,/,’ AVERAGE NUMBER OF’, 
$ » NEWTON-RAPHSON ITERATIONS’ ,T60,F7.3) 


C#NT THRU HERE 
END 


D Appendix Subroutine BESTEP 


c 
C 
C 
C 


aaagadnaaanagaagaagaagaagaagaanaandnanangngaagaagaanaaaa 


$ 
$ 


PAAAH OHA H 


NUMERICAL INTEGRATION OF 1 AND 2 DIMENSIONAL THETA 


SUBROUTINE ESTEP(QPT,COV,TDM,PMEAN ,QVAR,XMEANI,T, 


LK , LKTOT, GAMMA ,SIGMA, 
COVI ,GVEC ,GMU , GWN) 


IMPLICIT REAL*8 (A-H,0-Z) 
INCLUDE ’BIGMAT.FOR’ 
INCLUDE ’GLOBAL.FOR’ 


LOGICAL 
INTEGER*2 
REAL*4 
REAL*8 


QDIAG 
GVEC(NTOTG) 

T(NRANK) 

QPT(NQPT) , COV(NSCALE ,NSCALE) , TDM(NRANK , NSCALE) , 
PMEAN(NSCALE) ,QVAR(NSCALE) , 

XMEANI(NSCALE), — 

LK(NQPT, NSCALE) , LKTOT, 

GAMMA (NRANK , NSCALE) , SIGMA(NSCALE, NSCALE) , 
COVI(NSCALE,NSCALE) , 

GMU(NSCALE, TOTLEV) , 

GWN(TOTLEV) 


QPT GRID OF QUADRATURE POINTS CREATED IN PHASE1 


COV,TDM ACCUMALTES SUMS OF EXAMINEE POSTERIOR VARIANCES AND MEANS 


PMEAN 


QVAR 


HOLDS PRIOR MEANS COMPUTED FROM INP''TS OF PHASE1 WORK FILE 
AND GAMMA AND SIGMA FROM PREVIOUS ITERATION 


APPROXIMATE "VAR" OF EXAMINEE LIKELIHOOD FUNCTION 
INPUT FROM PHASE1 WORK FILE. USED ONLY IN DIAGNOSTIC 


REPORTING 


XMEANI 


LKTOT 


GAMMA , SIGMA 
PREVIOUS ITERATION 


EXAMINEE POSTERIOR MEAN 
HOLDS CONDITIONING VARIABLES READ FROM PHASE1 WORK FILE 


EXAMINEE LIKELIHOOD FUNCTION EVALUATED ON A GRID FOR 
EACH SCALE, INPUT FROM PHASE1 WORK FILE 


-2*L0G LIKELIHOOD FUNCTON 


VALUES OF GAMMA AND SIGMA FROM 


COVI EXAMINEE POSTERIOR VARIANCE 


GMU,GWN USED TO ACCUMULATE GROUP SUMMARY STATISTICS IF 
REQUESTED 
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aadnaggnganagagnagagaggagaagaggaagagangagagngnaagagaagagaagaagagaagaaa Qqaaa 


aaa 


9000 


VARIABLES USED TO CALCULATE POSTERIOR MEANS AND VARIANCES 


REAL*8 


LKCON 


RHO 


DET 


$1,S2 


P1,P2 


POSTI 


DIFF 


RHO, LKCON ,DET,B,S1,S2,P1,P2,POSTI , DIFF 


WIDTH OF QUADRATURE INTERVAL OR AREA OF SQUARE 
FOR NORMING LIKELIHOOD INTEGRAL CALCULATIONS 


FOR THE TWO DIMENSIONAL CASE, RHO IS THE PRIOR 
CORRELATION BETWEEN SUBSCALES 


DETERMINANT OF PRIOR VARIANCE SIGMA 

REGRESSION COEFFICIENT OF SCALE 2 ON SCALE 1 
USED FOR COMPUTING PRIOR DISTRIBUTION IN TwO 
STEPS 

CONSTANTS FOR COMPUTING PRIOR THAT 

INCLUDE NUMERICAL BOUND CHECKS 

PRIOR FOR SCALE 1, CONDITIONAL PRIOR FOR SCALE 2 
GIVEN SCALE 1 

POSTERIOR EVALUATED AT EACH QUADRATURE POINT 


NUMERATOR OF PRIOR 


THE VARIABLE LKTOT IS NOW ALWAYS UNIDIMENSIONAL 


CHARACTER*30 
CHARACTER*40 


DATA HEDCOV 
DATA HEDCTM 


SUBJID 

HEDCOV , HEDCTM 

/* ADJUSTED COVARIANCE MATRIX tf 
/*CHOLESKY TRANSFORMATION MATRIX / 


IF (NSCALE.GT.2) THEN 
WRITE (6, 9000) 
FORMAT(’?1 THIS ESTEP ONLY WORKS WITH 1 AND’ 


STOP 
END IF 


DSML = 1D-6 
DBIG = 14D1 


» 2 DIMENSIONAL THETA. ’) 


CONSTANTS NEEDED FOR POSTERIOR CALCULATIONS 
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LKCON=(QPMAX-QPMIN) /DFLOAT (NQPT-1) 
IF (NSCALE.EQ.2) LKCON=LKCON*LKCON 


IF (NSCALE.EQ.1)DET=SIGMA(1,1) 


IF (NSCALE.EQ.2) DET=SIGMA(1,1)*SIGMA(2,2)-SIGMA(1,2)*SIGMA(2,1) 


PUT IN NUMERIC BOUND AND CHECK FOR LESS THAN FULL RANK SIGMA 


ana 


IF (DET.LT.DSML) THEN 
DET=DSML 
CALL HEDOUT 
NLEFT=NLEFT-4 
WRITE(6,7100) ICYCL 
7100 FORMAT(’OTHE DETERMINANT OF SIGMA IS CLOSE TO ZERO. ’, 
$ » ICYCL=’,1I3,/,’ A NUMERIC BOUNDARY VALUE WAS ’, 
$ 7 INSERTED TO PREVENT PROGRAM TERMINATION.’,/, 
$ ’ OUTPUT THAT FOLLOWS IS UNRELIABLE. ’) 
END IF 


CONSTANTS FOR COMPUTING THE PRIOR DISTRIBUTION IN 
TWO STEPS 


Qaaaa 


S1=MAX( DSML,MIN(DBIG,2D0*SIGMA(1,1)) ) 

IF (NSCALE.EQ.2) THEN 
RHO=SIGMA(1,2) /DSQRT(SIGMA(1,1)*SIGMA(2,2)) 
B=RHO*DSQRT (SIGMA(2, 2) /SIGMA(1,1)) 

S2=MAX( DSML,MIN(DBIG, 2D0*(1-RHO*RHO)*SIGMA(2,2)) ) 

END IF 

Cc 
Cc INITIALIZE LKTOT, COV & TDM 
C*VDIR: ASSUME COUNT(2) 
LKTOT = ODO 
DO 40 ISCAL=1,NSCALE 
C*VDIR: ASSUME COUNT(2) 
DO 20 JSCAL=1 ,NSCALE 
COV(JSCAL,ISCAL) = ODO 
20 CONTINUE 
C*VDIR: ASSUME COUNT(50) 
DO 30 IRANK=1,NRANK 
TDMCIRANK,ISCAL) = ODO 
30 CONTINUE 
40 CONTINUE 
C*VDIR: ASSUME COUNT(10) 
DO 60 K=1,TOTLEV 
C*VDIR: ASSUME COUNT(2) 
DO 50 ISCAL=1 ,NSCALE 
GMUCISCAL,K) = ODO 


50 CONTINUE 
60 CONTINUE 
REWIND (IUNIT) 


ISUBJ = 0 
MAIN LOOP OVER SUBJECTS 


aaa 


C BESTEP 1. 
c 


100 READ(IUNIT,END=500) SUBJID(1:NIDW) ,KUSE,WGHT,MLO,MHI, 
$ (T(N) ,N=1,NRANK) , (QVAR(J) , J=1,NSCALE) , 
$ ((LK(IQPT,J) , IQPT=MLO,MHI) , J=1,NSCALE) , 
$ (GVEC(I) , l=1, NTOTG) 
IF (KUSE .EQ. 0) GO TO 100 


c 
C BESTEP 1. END 
c 
ISUBJ = ISUBJ+1 
QDIAG = ICYCL .EQ. 1 .AND. ISUBJ .LE. NQC 
C*VDIR: ASSUME COUNT(2) 
c 
c 
C BESTEP 2. 
c 


DO 120 ISCAL=1,NSCALE 
PMI = ODO 
C*VDIR: ASSUME COUNT(50) 
DO 110 IRANK=1,NRANK 
PMI = PMI+GAMMA(IRANK , ISCAL) *T(IRANK) 
110 CONTINUE 
PMEAN(ISCAL) = PMI 
120 CONTINUE 
Cc 


C BESTEP 2. END 


Cc 
IF (QDIAG) THEN 
CALL HEDOUT 
WRITE(6,8000) ISUBJ,ICYCL 
NLEFT = 55 
END IF 
LOOP OVER SCALES 
COMPUTE PRIOR & POSTERIOR DISTN’S 
ESTIMATES OF SUBJECT MEANS 
. & VARIANCES 
*VDIR: ASSUME COUNT(2) 


qgqaaqaqagngaa 


TOTAL=0D0 
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Q Q 7ea 


Q ? 


Q 


Q 


XMEANI (1) =0D0 

COVI(1,1)=0D0 

IF (NSCALE.EQ.2) THEN 
XMEANI(2)=0D0 
COVI(2,1)=0D0 
COVI(1,2)=0D0 
COVI(2,2)=0D0 

END IF 


COMPUTE PRIOR FOR SCALE 1 
BESTEP 3. 
DO 130 IQPT=MLO,MHI 


BESTEP 3(A) 


DIFF = (QPT(IQPT)-PMEAN(1))##2 
DIFF = MAX(DSML,MIN(DBIG,DIFF/S1) ) 
P1 = DEXP(-DIFF) 


BESTEP 3(A) END 


IF(NSCALE.EQ.2) THEN 


THE BIVARIATE PRIOR COMPUTED USING THE FACTORIZATION 
F (THETA1) *F (THETA2 | THETA1) 


DO 140 JQPT=MLO,MHI 
BESTEP 3(B) (TWO SCALES) 


DIFF = (QPT(JQPT)- (PMEAN(2)+B*(QPT(IQPT) 
$ -PMEAN(1))) )##2 

DIFF = MAX(DSML,MIN(DBIG,DIFF/S2)) 

P2=DEXP(-DIFF) 


BESTEP 3(C) (TWO SCALES) 
POSTI=P1*P2*LK (IQPT, 1) *LK (JQPT, 2) 
BESTEP 4. (TWO SCALES) 
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IMPLEMENTING EQUATIONS (4.1) AND (4.2) 


qgqaaa 


XMEANI (1) =XMEANI(1)+QPT(IQPT) *POSTI 

XMEANI (2) =XMEANI(2)+QPT(JQPT) *POSTI 
COVI(1,1)=COVI(1,1)+QPTCIQPT) *QPT(IQPT) *POSTI 
COVI(1,2)=COVI(1,2)+QPTCIQPT) *QPT(JQPT) *POSTI 
COVI(2,1)=COVI(1,2) 
COVI(2,2)=COVI(2,2)+QPT(JQPT) *QPT(JQPT) *POSTI 


Cc 
C BESTEP 5. (TWO SCALES) 
Cc 
TOTAL=TOTAL+POSTI 
140 CONTINUE 
Cc 
Cc 1 DIMENSION ONLY 
C 
ELSE 
Cc 
C BESTEP 3(C) (ONE SCALE) 
Cc 
POSTI=P1#*LK(IQPT, 1) 
Cc 
C BESTEP 4. (ONE SCALE) 
Cc 


XMEANI (1)=XMEANI (1)+QPT(IQPT) *POSTI 
COVI(1,1)=COVI(1,1)+QPT(IQPT) *QPT(IQPT) *POSTI 


C BESTEP 5. (ONE SCALE) 


TOTAL=TOTAL+POSTI 
END IF 
130 CONTINUE 
E 
C BESTEP 6. 
c 
C LIKELIHOOD CALCULATION: NOTE THAT LKCON AND DET 
C COMPUTED ONCE INITIALLY AND DIFFERENTLY FOR 1 AND 2 
C DIMENSIONAL CASES. COMPUTES -2* LOG OF LIKELIHOOD 
C FUNCTION IN (1.4) 
c 
c 
LKTOT = LKTOT+WGHT*( DLOG(TOTAL*LKCON)-0.5*DLOG(DET) ) 
c 
C BESTEP 7. 
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NORM AND CENTER POSTERIOR MOMENTS 


aaa 


XMEANI (1) =XMEANI (1) /TOTAL 

COVI(1,1)=COVI(1,1) /TOTAL 

COVI(1,1)=COVI(1,1)~XMEANI (1) *XMEANI (1) 

IF (NSCALE.EQ.2)THEN 
XMEANI (2) =XMEANI (2) /TOTAL 
COVI(2,2)=COVI(2,2)/TOTAL 
COVI(2,2)=COVI(2,2)-XMEANI (2) *XMEANI (2) 
COVI(2,1)=COVI(2,1)/TOTAL 
COVI(2,1)=COVI(2, 1)-XMEANI (2) *XMEANI (1) 
COVI(1,2)=COVI(2,1) 


END IF 
c 
C BESTEP 7. END 
c 
c 


IF (QDIAG) THEN 
WRITE(6 ,8030) 
NLEFT = NLEFT-3 
C*VDIR: ASSUME COUNT(2) 

DO 210 ISCAL=1,NSCALE 
QQVAR = ODO 
IF (QVAR(ISCAL) .GT. ODO) QQVAR = 1D0/QVARCISCAL) 
WRITE(6,8040) ISCAL,PMEANCISCAL) ,SIGMACISCAL, ISCAL) ,QQVAR, 

$ XMEANI CISCAL) , COVI(ISCAL, ISCAL) 
NLEFT = NLEFT-1 
210 CONTINUE 


END IF 
Cc 

IF (QDIAG) THEN 

CAL]. PRTMAT(COVI , NSCALE ,CM(KSNAM) ,NSCALE,CM(KSNAM) ,NSCALE, 

$ HEDCOV) 

END IF 
Cc 
Cc 
C BESTEP 8. 
Cc 
Cc ACCUMULATE TDM AND COV WITHOUT IMPUTATIONS ~ 
Cc EQUATIONS (1.7) AND (1.8) 
Cc 
Cc 

DO 340 JSCAL=1,NSCALE 
DO 330 IRANK=1,NRANK 
TDMCIRANK , JSCAL) = TDMCIRANK , JSCAL)+ 

$ WGHT*T CIRANK) *XMEANI (JSCAL) 

330 CONTINUE 
42 
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XMFAC=XMEANI ( JSCAL) -PMEAN( JSCAL) 


DO 350 ISCAL=JSCAL,NSCALE 
COVCISCAL,JSCAL) = COV(ISCAL, JSCAL)+ 


$ WGHT*COVI(ISCAL, JSCAL) 
COV(ISCAL,JSCAL) = COV(ISCAL, JSCAL)+ 
$ WGHT* (XMEANI (ISCAL) -PMEAN CISCAL) ) *XMFAC 
COV(JSCAL, ISCAL)=COV(ISCAL, JSCAL) 
350 CONTINUE 


340 CONTINUE 


c 
C BESTEP 8. END 


### REPLACE GROUP MEAN FROM IMPUTATION ESTIMATE WITH 
GROUP MEAN FROM ANALYTIC CALCULATION 
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DO 440 K=1,NTOTG 
ILEV = GVEC(K) 
IF (ILEV .GT. 0) THEN 
DO 430 ISCAL=1,NSCALE 
GMUCISCAL,ILEV) = GMUCISCAL, ILEV)+WGHT*XMEANI CISCAL) 

430 CONTINUE 

END IF 
440 CONTINUE 


GO TO 100 
500 CONTINUE 


RETURN -2LOGLIK AND NORM WEIGHTS TO ADD TO SAMPLE SIZE 


qaaaaa 


LKTOT=-2D0*LKTOT* (DFLOAT (ISUBJ) /WGTSUM) 


aQ 


DO 530 K=1,TOTLEV 
WT = GWN(K) 
IF (WI .NE. ODO) THEN 
DO 520 ISCAL=1,NSCALE 
GMUCISCAL,K) = GMUCISCAL,K)/WT 

520 CONTINUE 

END IF 
530 CONTINUE 


RETURN 
c 


Cc 
8000 FORMAT(’ <DIAG> PRINTOUT FOR SUBJECT NUMBER ’,14, 


43 


$ : CYCLE = ’,14) 
8010 FORMAT(’OSCALE’,13,’ THETA LK 
$ N.POST’) 
8020 FORMAT(’ ’,1I5,F10.3,4F10.5) 
8030 FORMAT(’?0 2 -=---= PRIOR--------- 
$ /? SCALE MU SIGMA VAR(LK) 
8040 FORMAT(’ ’,15,5F10.5) 
C 
C### NT FIXING ALTERNATE RETURN CALLS TO SYMINV 
Cc 
7000 WRITE(6,7010) 
7010 FORMAT(’ CALL TO SYMINV FAILED IN ESTEP’) 
STOP 
END 


PRIOR 


MEAN 


POST’, 


VAR ’) 


& 
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Appendix Subroutine NESTEP 


THIS ESTEP IS MODIFIED TO USE A NORMAL APPROXIMATION 
FOR MULTIVARIATE ESTEP. THE NORMAL APPROXIMATION IS 
DESCRIBED IN THE 1988 TECHNICAL REPORT (MGROUP AT THAT 
TIME DID NOT ACTUALLY IMPLEMENT THE NORMAL APPROXIMATION 
DESCRIBED THERE). ADDITIONAL DOCUMENTATION IS PROVIDED IN | 
“THE E-STEP OF THE EM ALGORITHM" 


SUBROUTINE ESTEP(COV,TDM,PMEAN,XMEANI ,QMEAN,QVAR,T, 

$ LKTOT, GAMMA, SIGMA, 

$ COVI ,SINV ,WMEAN , GVEC , GMU , GWN) 

IMPLICIT REAL*8 (A-H,0-Z) 

INCLUDE ’BIGMAT.FOR’ 

INCLUDE ’GLOBAL.FOR’ 

LOGICAL QDIAG 

INTEGER*2 GVEC(NTOTG) 

REAL*4 T(NRANK) 

REAL*8 COV(NSCALE,NSCALE) , TDM(NRANK ,NSCALE) , 
PMEAN(NSCALE) , XMEANI(NSCALE) , 
QMEAN (NSCALE) ,QVAR(NSCALE) , LKTOT, 
GAMMA (NRANK , NSCALE) , SIGMA(NSCALE,NSCALE) , 
COVI(NSCALE,NSCALE) ,SINV(NSCALE,NSCALE) , 
WMEAN (NSCALE) , GMU(NSCALE, TOTLEV) , 
GWN(TOTLEV) 


PARAAA SH 


VARIABLES FOR EVALUATING THE LIKELIHOOD FUNCTION 
REAL#8 DET , QFORM , QNORM , PI , DNSCAL 
CHARACTER*30 SUBJID 
CHARACTER*40 HEDCOV , HEDCTM 
PARAMETER (PI=3 . 141592653589793238) 


DATA HEDCOV /’ ADJUSTED COVARIANCE MATRIX >/ 
DATA HEDCTM /’CHOLESKY TRANSFORMATION MATRIX / 


QMEAN MEAN OF NORMAL APPROXIMATION TO EXAMINEE LIKELIHOOD 
READ FROM PHASE1 WORK FILE 

QVAR VAR OF NORMAL APPROXIMATION TO EXAMINEE LIKELIHOOD 
READ FROM PHASE1 WORK FILE. THE INVERSE OF THIS 
VARIANCE IS ACTUALLY READ FROM THE PHASE1 WORK FILE. 

PMEAN PRIOR MEAN COMPUTED WITH INPUT FROM PHASE1 WORK FILE 

GAMMA,SIGMA ESTIMATES. FROM THE PREVIOUS ITERATION 


XMEANI,COVI POSTERIOR MEAN AND VARIANCE FOR EACH EXAMINEE 


COVI IS REUSED TO HOLD A DIFFERENT MATRIX IN THE 
CALCULATION OF THE APPROXIMATE LIKELIHOOD FUNCTION. 


T CONDITIONING VARIABLES READ FROM PHASE1 WORK FILE 
GWN,GMU GROUP REPORTING VARIABLES 

LKTOT NEGATIVE TWICE THE LOG LIKELIHOOD 

SINV CONTAINS THE INVERSE OF SIGMA 

TDM,COV ACCUMULATE POSTERIOR MEANS AND VARIANCES 


WMEAN WORK SPACE FOR COMPUTING POSTERIOR MEANS 


NOTE THAT LKTOT IS NOW A SCALAR REGARDLESS OF 


Cc 
Cc 
Cc 
Cc 
Cc 
Cc 
Cc 
Cc 
Cc 
C 
Cc 
Cc 
Cc 
Cc 
Cc 
Cc 
C 
Cc 
Cc 
Cc 
C 


NSCALE 
DNSCAL=DFLOAT (NSCALE) 
DSML = 1D-6 
DBIG = 14D1 
Cc PUT SIGMA 
Cc IN SINV AND INVERT SINV 


C*VDIR: ASSUME COUNT (2) 
DO 10 JSCAL=1,NSCALE 
DO 15 ISCAL=1,NSCALE 
SINVCISCAL, JSCAL)=SIGMA(ISCAL, JSCAL) 


15 CONTINUE 
10 CONTINUE 
CALL SYMINV(SINV,NSCALE, DET, #7000) 
Cc 
Cc INITIALIZE LKTOT, COV & TDM 
C*VDIR: ASSUME COUNT(2) 


LKTOT=0D0 
DO 40 ISCAL=1,NSCALE 
C*VDIR: ASSUME COUNT(2) 
DO 20 JSCAL=1 ,NSCALE 
COV(JSCAL,ISCAL) = ODO 
20 CONTINUE 
C*VDIR: ASSUME COUNT(50) 
DO 30 IRANK=1,NRANK 
TDM(IRANK,ISCAL) = ODO 
30 CONTINUE 
40 CONTINUE 


Cc 
Cc INITIALIZE GROUP ACCUMULATION VARIABLES 
Cc 
Cc 


*VDIR: ASSUME COUNT(10) 
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LO 60 K=1, TOTLEV 
C*VDIR: ASSUME COUNT (2) 

DO 50 ISCAL=1,NSCALE 

GMUCISCAL,K) = ODO 
50 CONTINUE 
60 CONTINUE 
REWIND (IUNIT) 

ISUBJ = 0 


MAIN LOOP OVER SUBJECTS 


Q Qaaaa 


NESTEP 1. : 


READ STATEMENT MODIFIED TO READ ONLY INFO FROM 
APPROXIMATING NORMAL FOR THE LIKELIHOOD FROM 
READDF/READMF 


QaaQaaa 


100 READCIUNIT,END=500) SUBJID(1:NIDW) ,KUSE,WGHT, 
$ (T(N) ,N=1,NRANK) , (QMEAN(J) , J=1,NSCALE) , 
$ (QVAR(J) , J=1,NSCALE) , 
$ (GVEC(I) ,I=1,NTOTG) 


c ; 
IF (KUSE .EQ. 0) GO TO 100 
c 
C NESTEP 1. END 
c 
ISUBJ = ISUBJ+1 
c 
QDIAG = ICYCL .EQ. 1 .AND. ISUBJ .LE. NQC 
c 
C NESTEP 2. 
g 
c COMPUTE THE PRIOR MEAN FOR THETA 
c 


C*VDIR: ASSUME COUNT(2) 
DO 120 ISCAL=1,NSCALE 
PMI = ODO 
C*VDIR: ASSUME COUNT(50) 
DO 110 IRANK=1,NRANK 
PMI = PMI+GAMMA(IRANK , ISCAL) *T(IRANK) 
110 CONTINUE 
PMEANCISCAL) = PMI 
120 CONTINUE 
c ; 


C NESTEP 2. END 


Cc 
IF (QDIAG) THEN 
CALL HEDOUT 
WRITE(6,8000) ISUBJ,ICYCL 
NLEFT = 55 
END IF 
Cc 
C COMPUTE THE POSTERIOR VARIANCE OF THETA USING THE PRIOR VARIANCE 
Cc AND THE (DIAGONAL) VARIANCE OF THE NORMAL DISTRIBUTION 
C APPROXIMATING THE THETA-LIKELIHOOD. NOTE THAT QVAR 
Cc ALREADY CONTAINS THE INVERSE OF THE APPROXIMATING VARIANCE 
Cc 
Cc 
C NESTEP 3. 
Cc 
Cc 
Cc COMPUTE THE POSTERIOR VARIANCE OF THETA USING (4.5) 
Cc 


DO 130 JSCAL=1,NSCALL 
DO 140 ISCAL=1,NSCALE 
- COVICISCAL, JSCAL)=SINV(ISCAL, JSCAL) 
140 CONTINUE 
COVI CJSCAL, JSCAL)=COVI(JSCAL , JSCAL) +QVAR (JSCAL) 
130 CONTINUE 
Cc 
CALL SYMINV(COVI , NSCALE, DET, *7000) 


NESTEP 4. 


COMPUTE THE POSTERIOR MEAN OF THETA USING (4.6) 
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DO 150 JSCAL=1,NSCALE 
WMEAN ( JSCAL) =QVAR( JSCAL) *QMEAN ( JSCAL) 
DO 160 ISCAL=1 ,NSCALE 
; WMEAN ( JSCAL) =WMEAN (JSCAL) +SINV CISCAL, JSCAL) *PMEAN (ISCAL) 
160 CONTINUE 
150 CONTINUE 


DO 170 JSCAL=1,NSCALE 
XMEANI (JSCAL)=0D0 
DO 180 ISCAL=1,NSCALE 
XMEANI (JSCAL) =XMEANI (JSCAL) +COVI CISCAL, JSCAL) *WMEAN (ISCAL) 
180 CONTINUE 
170 CONTINUE 
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C NESTEP 4. END 
c 
g 
G 


IF (QDIAG) THEN 
WRITE(6 ,8030) 
NLEFT = NLEFT-3 
C*VDIR: ASSUME COUNT(2) 
DO 210 ISCAL=1,NSCALE 
QPVAR = ODO 
IF (QVAR(ISCAL) .GT. ODO) QPVAR = 1D0/QVAR(ISCAL) 
WRITE(6,8040) ISCAL,PMEAN(ISCAL) ,SIGMACISCAL, ISCAL) ,QPVAR, 
$ XMEANI(ISCAL) ,COVI(ISCAL, ISCAL) 
NLEFT = NLEFT-1 
210 CONTINUE 
END IF 


IF (QDIAG) THEN 

CALL PRTMAT(COVI ,NSCALE,CM(KSNAM) ,NSCALE, CM(KSNAM) ,NSCALE, 
$ HEDCOV) 
END IF 


c 
c 
C NESTEP 5. 


ACCUMULATE TDM AND COV WITHOUT IMPUTATIONS 
EQUATIONS (1.7) AND (1.8) 
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DO 340 JSCAL=1,NSCALE 
DO 330 IRANK=1 ,NRANK 
TDMCIRANK, JSCAL) = TDMCIRANK, JSCAL)+ 
$ WGHT*T (IRANK) *XMEANI ( JSCAL) 
330 CONTINUE 
XMFAC=XMEANI ( JSCAL) -PMEAN (JSCAL) 


DO 350 ISCAL=JSCAL ,NSCALE 
COVCISCAL, JSCAL) = COV(ISCAL, JSCAL)+ 
$ WGHT*COVI (ISCAL, JSCAL) 
COV(ISCAL, JSCAL) = COV(ISCAL, JSCAL)+ 
$ WGHT* (XMEANI (ISCAL) -PMEAN(ISCAL) ) *XMFAC 
COV(JSCAL, ISCAL)=COVCISCAL, JSCAL) 
350 CONTINUE 
340 CONTINUE 


##8 END OF CODE TO COMPUTE TDM AND COV 


ana -aAA 


49 


QU 


NESTEP 6. 


C 
c 
Cc ###CODE TO COMPUTE THE APPROXIMATE (WEIGHTED) LOG LIKELIHOOD 
Cc COVI IS REUSED TO HOLD A VARIANCE MATRIX 

c 

QFORM=0D0 

QNORM=(DNSCAL/2D0) *DLOG (2*PI) 

IF (NSCALE.EQ.1)THEN 
DET=SIGMA(1,1)+(1.0/MAX (DSML ,QVAR(1))) 
COVI(1,1)=1.0/DET 
QNORM=QNORM+.5*DLOG(1.0/MAX(DSML , QVAR(1))) 

QFORM=( QMEAN(1)-PMEAN(1) )*( QMEAN(1)- 
$ PMEAN(1) )*COVI(1,1) 
ELSE 
DO 360 JSCAL=1,NSCALE 
QNORM=QNORM+ .5*DLOG(1.0/MAX(DSML ,QVAR(JSCAL) )) 
DO 370 ISCAL=1,NSCALE 
COVICISCAL, JSCAL) =SIGMA(ISCAL, JSCAL) 
370 CONTINUE 
COVI(JSCAL, JSCAL)=COVI (JSCAL, JSCAL) + 
$ (1.0/MAX(DSML, QVAR(JSCAL) )) 
360 CONTINUE 


CALL SYMINV(COVI ,NSCALE , DET, +7000) 
DET=DEXP (DET) 


DO 380 JSCAL#1, (NSCALE-1) 
QFORM=QFORM+ 
$ ( QMEANCJSCAL)-PMEANCJSCAL) )*( QMEAN(JSCAL)- 
$ PMEAN(JSCAL) )*COVI(JSCAL, JSCAL) 
DO 390 ISCAL=(JSCAL+1) ,NSCALE 
QFORM=QFORM+ 
$ 2.0*( QMEANCISCAL)-PMEAN(ISCAL) )*( QMEAN(JSCAL)- 
$ PMEAN(JSCAL) )*COVI(ISCAL, JSCAL) 
390 CONTINUE 
380 CONTINUE 
QFORM=QFORM+( QMEAN(NSCALE)-PMEAN(NSCALE) )*( QMEAN(NSCALE) - 
$ PMEAN(NSCALE) )*COVI(NSCALE,NSCALE) 
END IF 


Q 


LKTOT=LKTOT+WGHT*( QNORM-0.5*(LOG(DET)+QFORM) ) 


c 
C NESTEP 6. END 


##% REPLACE GROUP MEAN FROM IMPUTATION ESTIMATE WITH 
GROUP MEAN FROM ANALYTIC CALCULATION 


Qaaaaa 


DO 440 K=1,NTOTG 


06 


ILEV = GVEC(K) 
IF (ILEV .GT. 0) THEN 
DO 430 ISCAL=1,NSCALE 
GMUCISCAL, ILEV) = GMUCISCAL, ILEV)+WGHT*XMEANI (ISCAL) 
430 CONTINUE 
END IF 
440 CONTINUE 


GO TO 100 
500 CONTINUE 


C 
Cc 
Cc RETURN -2LOGLIK AND NORM WEIGHTS TO ADD TO SAMPLE SIZE 
Cc 
Cc 


LKTOT=-2D0*LKTOT* (DFLOAT (ISUBJ) /WGTSUM) 


DO 530 K=1,TOTLEV 
WI = GWN(K) 
IF (WT .NE. ODO) THEN 
DO 520 ISCAL=1,NSCALE 
GMUCISCAL,K) = GMUCISCAL,K)/WT 

520 CONTINUE 

END IF 
530 CONTINUE 


RETURN 
Cc 
Cc 
8000 FORMAT(’ <DIAG> PRINTOUT FOR SUBJECT NUMBER ’,14, 
$ 4 CYCLE = ’,14) 
8010 FORMAT(’OSCALE’,I3,’ THETA LK PRIOR POST’, 
$ N.POST’) 
8020 FORMAT(’ ’,15,F10.3,4F10.5) 
86:0 FORMAT(’?0 2 -=-=- PRIOR<---nn---- tee POSTERIOR-----’ 
$ /’ SCALE MU SIGMA VAR(LK) MEAN VAR ’) 
8040 FORMAT(’ ’,1I5,5F10.5) 
8060 FORMAT(’OIMPUTATION #’ ,14) 
8070 FORMAT(’ Z: ’,10F12.4) 
8080 FORMAT(’ D: ’,10F12.4) 
Cc 
C### NT FIXING ALTERNATE RETURN CALLS TO SYMINV 
c 
7000 WRITE(6,7010) 
7010 FORMAT(’ CALL TO SYMINV FAILED IN ESTEP’) 
STOP 
END 


F Appendix Subroutine CESTEP 


Cc 
Cc 
Cc 
Cc 
Cc 


aaa qaaaaa 
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THIS ESTEP IS MODIFIED TO USE HIGHER ORDER ASYPMTOTIC 
CORRECTIONS TO APPROXIMATE EXAMINEE POSTERIOR MEANS AND 
VARIANCES 


SUBROUTINE ESTEP(COV,TDM,PMEAN,XMEANI ,QMEAN,QVAR,T, 
LKTOT , GAMMA ,SIGMA, XSTEP, 
COVI ,SINV, WMEAN ,D,GVEC,GMU ,GWN,SLOPE, 
THRESH , GUESSC , CATEG , LINK ,NCAT , XITEM, 
DERIV,G,XBACK ,COVIB ,XMODEI , ISUBJ, 
IFAIL, INOTPD, IFCORM, IFCORV , AVEIT) 
IMPLICIT REAL*8 (A-H,0-Z) 
INCLUDE ’BIGMAT.FOR’ 
INCLUDE ’GLOBAL.FOR’ 
LOGICAL QDIAG 
INTEGER*2 GVEC(NTOTG) 
REAL*4 TCNRANK) 
REAL*8 COV(NSCALE , NSCALE) , TDM(NRANK , NSCALE) , 
PMEAN (NSCALE) , XMEANI (NSCALE) , XBACK (NSCALE) , 
QMEAN (NSCALE) ,QVAR(NSCALE) , 
LKTOT, 
GAMMA (NRANK , NSCALE) , SIGMA (NSCALE , NSCALE) , 
XSTEP(NSCALE) , 
COVI(NSCALE,NSCALE) , SINV(NSCALE,NSCALE) , 
WMEAN (NSCALE) ,D(NSCALE) ,GMU(NSCALE, TOTLEV) , 
GWN(TOTLEV) , DERIV(NSCALE,4) ,G(NSCALE,NSCALE) , 
COVIB(NSCALE , NSCALE) , XMODEI (NSCALE) 


PAAA SH 


PAAAAAGAGSHA 


ITEM PARAMETERS AND RESPONSES PASSED TO CESTEP SO DERIVATIVES 
OF THE LIKELIHOOD FUNCTION CAN BE COMPUTED DIRECTLY 


REAL*8 SLOPE(NITM) , THRESH(NITM) ,GUESSC(NITM) , CATEG (MAXCAT, NITM) 
REAL*4 XITEM(NITM) 
INTEGER*4 LINK(NITM) ,NCAT(NITM) 
VARIABLES FOR EVALUATING THE LIKELIHOOD FUNCTION 
REAL*8 DET ,QFORM, PI ,D2NSCA ,DLSIG ,DETG , LKNORM 
LKNORM COMPUTED IN READDF/READMF AND PASSED TO CESTEP. 
MAXIMUM VALUE OF THE EXAMINEE LOG LIKELIHOOD USED 
TO NORM THE LOG LIKELIHOOD FUNCTION SO THAT IT IS 


ROUGHLY COMPARABLE TO THE OTHER METHODS OF 
APPROXIMATING THE LOG LIKELIHOOD FUNCTION. 


CHARACTER*30 SUBJID 
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CHARACTER*40 HEDCOV , HEDCTM 

PARAMETER (PI=3. 141592653589793238) 

DATA HEDCOV /’ ADJUSTED COVARIANCE MATRIX <3 
DATA HEDCTM /*CHOLESKY TRANSFORMATION MATRIX f 


FOR THE NORMAL APPROXIMATION THAT IS CALCULATED FIRST 

QMEAN(J) AND QVAR(J) HOLD THE MEAN AND (INVERSE OF THE) VARIANCE OF 
THE APPROXIMATING NORMAL DISTRIBUTION OF THE EXAMINEE LIKELIHOOD 
FUNCTION FOR SCALE J. 

PMEAN HOLDS THE PRIOR MEANS 

SIGMA HOLDS THE PRIOR VARIANCES (THERE IS NO PVAR). 

XMEANI,COVI HOLD THE POSTERIOR MEAN AND VARIANCE. 

COVI IS REUSED TO HOLD A DIFFERENT MATRIX IN THE 

CALCULATION OF THE APPROXIMATE LIKELIHOOD FUNCTION. 


T HOLDS THE 
CONDITIONING VARIABLES READ FROM THE PHASE1 WORK FILE. 


TDT AND COV 
ACCUMULATE THE SUMS OF THE MEANS AND VARIANCES. 


XSTEP IS MAXIMUM 


CHANGE ALLOWED IN NEWTON-RAPHSON ITERATIONS AND MAXIMUM CHANGE 
ALLOWED IN CORRECTED APPROXIMATIONS. 


“ 
XMODEI IS THE POSTERIOR MODE. 
G INVERSE OF SECOND PARTIALS. 


XBACK,COVIB BACKUP OF NGROUP 
ESTIMATES USED IN CASE CORRECTIONS FAIL. 


SINV INVERSE OF SIGMA. 
D IS WORK VECTOR FOR NEWTON-RAPHSON ITERATIONS. 


NOTE THAT LKTOT IS NOW A SCALAR REGARDLESS OF 
NSCALE 


VARIABLES USED TO FIND MODAL VALUE OF THE POSTERIOR DISTRIBUTION 
BY NEWTON-RAPHSON 


REAL*8 XIMIN ,XIMAX, AVEIT , XLNLIK 
INTEGER IPD, INOTPD, IFAIL, IFCORM, IFCORV,ICFM,ICFV 


FROM COMMON XNR1,XNR2,XNR3,XNR4,XNRS,XNR6,MAXIT 


og 


Cc 
c 
Cc 
Cc 
c 
Cc 
c 
Cc 
Cc 
C 
Cc 
C 
c 
Cc 
C 
Cc 
C 
Cc 
Cc 
Cc 
c 
c 
c 
c 
c 
Cc 
Cc 
Cc 
Cc 
Cc 
c 
c 
c 
c 
Cc 
Cc 
c 
c 
Cc 
Cc 
Cc 
c 
Cc 
Cc 
c 
c 
c 
Cc 
c 
c 
c 


XNR1 PROPORTION OF CHANGE (MEASURED BY SD) FROM NORMAL 
APPROXIMATION TO THE LIKELIHOOD ALLOWED AT EACH STEP OF N-R 


XNR2 STOPPING CRITERIA APPLIED TO EACH COORDINATE IN THE 
NEWTON-RAPHSON ROUTINE 


XNR3 MAGNITUDE OF THE RATIO OF THE INFO MATRIX OBTAINED AT EACH 
N-R STEP AND THE CORRESPONDING APPROXIMATE INFORMATION FROM 
THE NORMAL APPROXIMATION TO THE LIKELIHOOD FUNCTION. 
THIS RATIO IS USED TO INSURE A POSITIVE DEFINITE MATRIX 
IN THE NEWTON RAPHSON ROUTINE 

XNR4 NUMBER OF XSTEPS OF TOTAL CHANGE ALLOWED IN XMEANI 


XNRS MAXIMUM RATIO OF MODAL VARIANCE ESTIMATE TO NORMAL 
BASED ESTIMATE ALLOWED 


XNR6 MINIMUM RATIO OF CORRECTED VARIANCES TO NORMAL BASED VARIANCES 


XIMIN,XIMAX MINIMUM AND MAXIMUM DIAGONAL ELEMENTS ALONG DIAGONAL 
OF THE NORMAL BASED INFORMATION MATRIX 


XLNLIK EXAMAINEE LOG LIKELIHOOD FUNCTION EVALUATED AT 
THE CONVERGED VALUE OF THE NEWTON ROUTINE FOR USE 
IN COMPUTING THE LIKELIHOOD. FUNCTION 

MAXIT MAXIMUM NUMBER OF ITERATIONS ALLOWED 


DERIV HOLDS THE FIRST 4 DERIVATIVES OF THE 
NEGATIVE LOGLIKELIHOOD FUNCTION 


G HOLDS SECOND DERIVATIVES OF THE NEGATIVE LOG POSTERIOR 


INOTPD COUNTS THE NUMBER OF TIMES THE INFO MATRIX IS NOT WELL 
IPD CONDITIONED 


XNLIK THETA LOG LIKELIHOOD 
AVEIT AVERAGE NUMBER OF ITERATIONS FOR THE NEWTON SEARCH 
IFAIL COUNTS THE NUMBER OF NEWTON SEARCH FAILURES 


IFCORM COUNTS THE NUMBER OF TIMES NO CORRECTION TO THE MEAN WAS 
ICFM APPLIED(WHEN THE NEWTON WAS SUCCESSFUL) 


IFCORV COUNTS THE NUMBER OF TIMES NO CORRECTION 


ICFV TO THE VARIANCE WAS APPLIED(WHEN THE 
NEWTON WAS SUCCESSFUL) 


60 


IFAIL=0 
IFCORM=0 
IFCORV=0 
INOTPD=0 
AVEIT=0D0 


D2NSCA=DFLOAT (NSCALE) /(2.0D0) 
DSML = 1D-6 
DBIG = 14D1i 


PUT SIGMA 
IN SINV AND INVERT SINV 


aqaaaa 


*VDIR: ASSUME COUNT(2) 
DO 10 JSCAL=1,NSCALE 
DO 15 ISCAL=1,NSCALE 
SINVCISCAL, JSCAL) =SIGMA(ISCAL, JSCAL) 


15 CONTINUE 
10 CONTINUE 
CALL SYMINV(SINV,NSCALE, DLSIG , +7000) 
c 
c INITIALIZE LKTOT, COV & TDM 
C*VDIR: ASSUME COUNT (2) 


LKTOT=0D0 
DO 40 ISCAL=1,NSCALE 
C*VDIR: ASSUME COUNT (2) 
DO 20 JSCAL=1,NSCALE 
COV(JSCAL,ISCAL) = ODO 
20 CONTINUE 
C*VDIR: ASSUME COUNT(50) 
DO 30 IRANK=1 ,NRANK 
TDMC(IRANK,ISCAL) = ODO 
30 CONTINUE 
40 CONTINUE 


Cc 
c INITIALIZE GROUP ACCUMULATION VARIABLES 
Cc 
Cc 


*VDIR: ASSUME COUNT(10) 
DO 60 K=1,TOTLEV 
C*VDIR: ASSUME COUNT (2) 
DO SO ISCAL=1 ,NSCALE 
GMUCISCAL,K) = ODO 


50 CONTINUE 
60 CONTINUE 
REWIND (IUNIT) 
ISUBJ = 0 
c 
c MAIN LOOP OVER SUBJECTS 
c 
C CESTEP 1. 


55 


READ STATEMENT MODIFIED TO READ ONLY INFO FROM 
APPROXIMATING NORMAL FOR THE LIKELIHOOD FROM 
READDF/READMF. IN ADDITION, THE CORRECTED NORMAL 
READS THE ITEM RESPONSES FOR COMPUTATION OF DERIVATIVES 
OF THE LOG LIKELIHOOD FUNCTION 


qgqaaaaaaa 


100 READ(IUNIT,END=500) SUBJID(1:NIDW) ,KUSE,WGHT, 
$ (T(N) ,N=1,NRANK) , (QMEAN(J) , J=1,NSCALE) , 
$ (QVAR(J) , J=1,NSCALE) , LKNORM, 
$ (GVEC(I) , I=1,NTOTG) , (XITEM(ITM) , ITM=1,NITM) 


c 
IF (KUSE .EQ. 0) GO TO 100 
c 
C CESTEP 1. END 
G 
ISUBJ = ISUBJ+1 
c 
QDIAG = ICYCL .EQ. 1 .AND. ISUBJ .LE. NQC 
c 
c 
c 
C CESTEP 2. 
c 
c COMPUTE THE PRIOR MEAN FOR THETA 
c 
C*VDIR: ASSUME COUNT(2) 


DO 120 ISCAL=1,NSCALE 
PMI = ODO 
C*VDIR: ASSUME COUNT(50) 
DO 110 IRANK=1,NRANK 
PMI = PMI+GAMMA (IRANK , ISCAL) *T(IRANK) 
110 CONTINUE 
PMEANCISCAL) = PMI 
120 CONTINUE 
c 


C CESTEP 2. END 


Cc 
IF (QDIAG) THEN 
CALL HEDOUT 
WRITE(6,8000) ISUBJ,ICYCL 
NLEFT = 55 
END IF 


qaaaQaQ 
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C 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


140 


130 


aaaaaa 


2010 


aa a 


aaa 


180 


CESTEP 3. 


USE NORMAL APPROXIMATION FROM 1990 TECHNICAL REPORT 
TO OBTAIN INITIAL STARTING VALVE TO OBTAIN THE MODE 
OF THETA 


COMPUTE THE POSTERIOR VARIANCE OF THETA USING THE PRIOR VARIANCE 
AND THE (DIAGONAL) VARIANCE OF THE NORMAL DISTRIBUTION 
APPROXIMATING THE EXAMINEE LIKELIHOOD. NOTE THAT QVAR 

ALREADY CONTAINS THE INVERSE OF THE APPROXIMATING VARIANCE 


DO 130 JSCAL=1,NSCALE 
DO 140 ISCAL=1,NSCALE 
COVICISCAL, JSCAL) =SINV(ISCAL, JSCAL) 
CONTINUE 
COVI(JSCAL, JSCAL)=COVI( JSCAL, JSCAL) +QVAR(JSCAL) 
CONTINUE 


COMPUTE THE MAX AND MIN DIAGONAL ELEMENTS. THESE BOUNDS ARE 
USED BY THE CODE TO CORRECT THE NORMAL APPROXIMATION. 


XIMIN=COVI(1,1) 
XIMAX=XIMIN 
DO 2010 JSCAL=2,NSCALE 
IF( COVICJSCAL, JSCAL) .GT.XIMAX )XIMAX=COVI(JSCAL, JSCAL) 
IF( COVICJSCAL, JSCAL) .LT.XIMIN )XIMIN=COVI(JSCAL, JSCAL) 
CONTINUE 


CALL SYMINV(COVI ,NSCALE, DET, #7000) 
COMPUTE THE POSTERIOR MEAN OF THETA 


DO 150 JSCAL=1,NSCALE 
WMEAN ( JSCAL) =QVAR( JSCAL) *QMEAN (JSCAL) 
DO 160 ISCAL=1,NSCALE 
WMEAN ( JSCAL) =WMEAN ( JSCAL) +SINV(ISCAL, JSCAL) *PMEAN (ISCAL) 
CONTINUE 
CONTINUE 


DO 170 JSCAL=1,NSCALE 
XMEANI (JSCAL) =0D0 
DO 180 ISCAL=1,NSCALE 
XMEANI (JSCAL) =XMEANI ( JSCAL) +COVI CISCAL, JSCAL) *WMEAN (ISCAL) 
CONTINUE 
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170 CONTINUE 
Cc 


C CESTEP 3. END 


Cc 
Cc 


CeeeeeeeeeeeSSKEKSEEEEERESESEREKER ERASER ERE REESE REREESEEEEEEERE SESE 


C NEWTON-RAPHSON ROUTINE TO FIND THE MODAL VALUE OF THE POSTERIOR 


WMEAN REUSED AS A WORK VECTOR 


CESTEP 4. 


COMPUTE THE MAX STEP SIZES ALLOWED FOR EACH COORDINATE. 
MAKE BACKUP COPY OF POSTERIOR MEANS AND VARIANCES FROM THE NORMAL 
APPROXIMATION IN CASE THE CORRECTION FAILS. 


THE STARTING VALUES OF NEGATIVE MEANS ARE INCREASED 
BECAUSE THE DISTRIBUTIONS ARE TYPICALLY SKEWED WITH 
THE MEAN MORE NEGATIVE THAN THE MODE 


aagaaaaaaaa Q eaaa 


DO 2020 JSCAL=1,NSCALE 

XSTEP ( JSCAL) =XNR1*DSQRT(COVI(JSCAL, JSCAL) ) 

XBACK ( JSCAL) =XMEANI (JSCAL) 

XMODEI (JSCAL) =XMEANI (JSCAL) 

IF (XMODEI (JSCAL) .LT.ODO) 

$ XMODEI ( JSCAL) =XMODEI ( JSCAL) +0 . SD0*XSTEP ( JSCAL) 

COVIB(JSCAL ,, JSCAL) =COVI (JSCAL, JSCAL) 

DO 2025 ISCAL=(JSCAL+1) ,NSCALE 
COVIBCISCAL, JSCAL)=COVI CISCAL , JSCAL) 
COVIB(JSCAL, ISCAL) =COVI (JSCAL, ISCAL) 

2025 CONTINUE 
2020 CONTINUE 


c 
C CESTEP 4. END 


c 
c 
c 
C SET CONSTANTS FOR ENSURING A WELL CONDITIONED INFO MATRIX 
c 
B2=XNR3*XIMAX 
DMIN=XIMIN/XNR3 
c 
C CESTEP 5. 
c 
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Cc BEGIN N-R ITERATIONS 
Cc 

IPD=0 

ITER=0 


2100 ITER=ITER+1 


COMPUTE THE DERIVATIVES OF THE LOG LIKELIHOOD FUNCTION 


Cc 
C CESTEP 5(A) 
Cc 
Cc 
Cc 


CALL DERVS2(XITEM,NITM,NSCALE,MAXCAT,SLOPE, THRESH, 
$ GUESSC , CATEG , LINK ,NCAT, XMODEI , DERIV , DFAC) 


COMPUTE THE DERIVATIVES OF THE NEGATIVE LOG POSTERIOR 


aaa 


DO 2210 JSCAL=1,NSCALE 
DERIV(JSCAL, 1)=DERIV(JSCAL, 1) 
DO 2220 ISCAL=1,NSCALE 
GCISCAL, JSCAL)=SINVCISCAL, JSCAL) 
DERIV(JSCAL, 1)=DERIV(JSCAL, 1)+ 
$ SINVCISCAL, JSCAL)*( XMODEICISCAL)-PMEANCISCAL) ) 
2220 CONTINUE 
G(JSCAL, JSCAL)=G(JSCAL, JSCAL) +DERIV(JSCAL, 2) 
2210 CONTINUE 


Cc 
C CESTEP 5(B) 
Cc 
Cc COMPUTE THE CHOLESKY DECOMPOSITION OF G 
Cc 
CALL DCHOL(G,NSCALE, NSCALE, INFO,B2,DMIN) 
IPD=IPD+INFO 
Cc 
Cc COMPUTE G*(-1)%*(FIRST DERIVATIVES) 
Cc 
CALL CHOLMU(G,DERIV(1,1) ,NSCALE, WMEAN ,D) 
Cc 
C CESTEP 5(C) 
Cc 
Cc UPDATE THE ESTIMATE OF ‘THE POSTERIOR MODE AND CHECK FOR 
Cc CONVERGENCE 
Cc 
ICONV=1 


DO 2230 JSCAL=1,NSCALE 
IF( DABS(D(JSCAL)) .GT.XSTEP(JSCAL) ) 
$ D(JSCAL)=SIGN( XSTEP(JSCAL) ,D(JSCAL) ) 
IF( DABS(D(JSCAL)) .GT.XNR2 ) ICONV=0 
XMODEI ( JSCAL) =XMODEI ( JSCAL) -D(JSCAL) 


2230 
Cc 


CONTINUE 
IF( (ITER.LT.MAXIT) .AND. CICONV.EQ.0) ) GO TO 2100 
COUNT AND REPORT THE AVERAGE NUMBER OF ITERATIONS REQUIRED 
AVEIT=AVEIT+ITER 


COUNT THE NUMBER OF SUBJECTS THAT HAD A NON POSITIVE 
DEFINITE MATRIX DURING SOME ITERATION 


IF (IPD.GT.0) INOTPD=INOTPD+1 


END OF N-R ITERATION 


CHSSEKAREKSEESEESSESESEEASEKSKSE SHEARS SSSA SSKAEKEREA EKER SEKEEERE EE 


Cc 


C CESTEP 5. END 


Cc 


CHHARERKERRESEREEREEEERAEESAEKESAAESEAKESRSE RAR REAEAEREESEERERERESE RE 


Cc 


C CESTEP 6. 


qaagqaaaa 


Q Q e@eaeanaa 


qaaa 


aaa 


APPLY HIGHER ORDER ASYMPTOTIC CORRECTIONS. 

IF THE N-R CODE FAILED TO ACHIEVE CONVERGENCE, THEN DO 
NOT PERFORM ADJUSTMENT, USE THE ORGINAL VALUES FROM THE 
NORMAL APPPROXIMATION TO THE LIKELIHOOD FUNCTION. 


ICFM=0 
ICFV=0 
IF CICONV.EQ. 1) THEN 


COMPUTE THE INVERSE OF THE FINAL INFORMATION MATRIX FROM 
THE N-R ROUTINE. THE VARIANCES THIS YIELDS ARE USED IN THE 
HIGHER ORDER CORRECTION TERMS. THE DETERMINANT OF THE 
ASYMPTOTIC VARIANCE MATRIX IS COMPUTED FOR USE IN THE 
CALCULATION OF THE LIKELIHOOD FUNCTION 


CESTEP 6(A) 


CALL CHOLIN(G,NSCALE,DETG) 

CHOLIN RETURNS THE DETERMINANT OF G BEFORE INVERSION. 
DETG=1D0/DETG 

FIRST COMPUTE THE HIGHER (3 AND 4) 
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Cc 


DERIVATIVES OF THE -LOG POSTERIOR(=DERVS OF -LOG LIK) AT 
THE FINAL ESTIMATE OF THE MODE FROM THE N-R. 

THE FIRST TWO DERIVATIVES ARE REUSED FROM THE 

N-R CALCULATIONS. THEY ARE NOT RECOMPUTED HERE. 


CALL DERVS(XITEM,NITM,NSCALE,MAXCAT ,SLOPE, THRESH ,GUESSC, 


CATEG , LINK , NCAT ,XMODEI , DERIV , DFAC,, XLNLIK) 


C CESTEP 6(B) 


C 


Cc 
Cc 
Cc 


aaa 


C 
C 
C 
c 
c 
c 
c 
c 


DO 2350 JSCAL=1,NSCALE 


FIRST ORDER CORRECTION TO THE MEANS IN (3.3) 
W=0D0 
DO 2375 ISCAL=1,NSCALE 
W=W+G(ISCAL, ISCAL) *G(ISCAL, JSCAL) *DERIV(ISCAL, 3) 


CONTINUE 
W=0.5D0*W 


XMEANI ( JSCAL) =XMODEI (JSCAL) -W 


CESTEP 6(C) 


CHECK SIZE OF THE CORRECTION TO THE MEAN 
BOUND AMOUNT OF CHANGE ALLOWED 
XS AND XD ARE LOCAL WORK VARIABLES 


XS=XNR4*XSTEP (JSCAL) 
XD=XMEANI ( JSCAL) -XBACK ( JSCAL) 


IF( DABS(XD) .GT. XS)THEN 


ICFM=1 
XMEANTI ( JSCAL)=XBACK ( JSCAL) +SIGN (XS , XD) 
END IF 


CESTEP 6(D) 


AVERAGE THE TWO VARIANCE ESTIMATES FOR 

BOUND CHECK. THE WEIGHTING OF THE AVERAGE 

WAS DETERMINED EMPIRICALLY FROM A NUMBER OF EXAMPLES. 
IF THE MODAL VARIANCE ESTIMATE IS TOO LARGE, THEN 
JUST USE THE NORMAL-BASED ESTIMATE 


IF( GCJSCAL, JSCAL) /COVIB(JSCAL,JSCAL) .GT. XNRS)THEN 
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W=COVIB(JSCAL, JSCAL) 


ELSE 
W=0..25D0*COVIB(JSCAL, JSCAL) +0. 75D0*G(JSCAL, JSCAL) . 
END IF 
Cc 
DO 2385 ISCAL=JSCAL,NSCALE . 


Cc 
C CESTEP 6(E) 


APPLY SECOND ORDER CORRECTION TO THE VARIANCES IN (3.4) 


aaaa 


W3=0D0 

W3A=0D0 

W4=0D0 

DO 2410 K=1,NSCALE 

W4=W4+DERIV(K,4)*G(K, ISCAL) *G(K , JSCAL) *G(K,K) 
W3A=W3A+DERIV(K ,3) *DERIV(K ,3) *G(K, ISCAL) * 
$ G(K, JSCAL) *G (K ,K) *G(K ,K) 
W3B=0D0 
DO 2420 L=1, (K-1) 
W3B=W3B+DERIV(L,3)*G(L,K)* 
( G(L, ISCAL)*G(L,JSCAL)*G(K,K) + 
G(L, ISCAL) *G(K, JSCAL) *G(L,K) + 
G(K, ISCAL) *G(L, JSCAL)*G(L,K) + 
G(K, ISCAL) *G(K, JSCAL)*G(L,L) ) 
2420 CONTINUE 
W3=W3+W3B*DERIV(K,3) 
2410 CONTINUE 
COVICISCAL, JSCAL)=GCISCAL, JSCAL) -0.SD0*W4+W3A+0. 5D0*W3 
COVI(JSCAL, ISCAL)=COVI CISCAL, JSCAL) 

2385 CONTINUE 


PAPAAA 


BOUND THE MAGNITUDE OF THE CORRECTION 


Cc 
C CESTEP 6(F) 
Cc 
Cc 
Cc 


IF(COVI(JSCAL,JSCAL)/W .GT. XNRS) THEN 


ICFV=1 
COVI(JSCAL, JSCAL) =W*XNRS 
ELSE IF(COVI(JSCAL, JSCAL) .LE.W*XNR6) THEN 
ICFV=1 
COVI(JSCAL, JSCAL)=W*XNR6 
END IF 
Cc 
c . 
2350 CONTINUE 
Cc 
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Cc USE CORRELATIONS FROM THE ORIGINAL NORMAL APPROXIMATION 
c IF THE CORRECTION FAILED ON ANY VARIANCE 
c 


IFCICFM.EQ. 1) IFCORM=IFCORM+1 
c 


C CESTEP 6(G) 


C 
IFCICFV.EQ.1) THEN 
IFCORV=IFCORV+1 
DO 2310 JSCAL=1 ,NSCALE 
DO 2305 ISCAL=(JSCAL+1) ,NSCALE 
XSD=DSQRT(COVI (JSCAL,, JSCAL) *COVI CISCAL, ISCAL) ) 
XBSD=DSQRT (COVIB( JSCAL, JSCAL) *COVIB(ISCAL , ISCAL) ) 
XBCOR=COVIB(ISCAL , JSCAL) /XBSD 
COVI CISCAL, JSCAL) =XSD*XBCOR 
COVI (JSCAL, ISCAL) =COVI (ISCAL, JSCAL) 
2305 CONTINUE 
2310 CONTINUE 
END IF 


CESTEP 6(G) END 


NEWTON-RAPHSHON FAILED TO CONVERGE, USE NORMAL-BASED MOMENTS 
ALREADY STORED IN XMEANI AND COVI 


aAaAaAaAaA-AnRAaAA Q ? 


ELSE 
IFAIL=IFAIL+1 
END IF 


aa aA 


IF (QDIAG) THEN 
WRITE(6 , 8030) 
NLEFT = NLEFT-3 
C*VDIR: ASSUME COUNT(2) 
DO 210 ISCAL#1,NSCALE 
QPVAR = ODO 
IF (QVARC(ISCAL) .GT. ODO) QPVAR = 1D0/QVAR(ISCAL) 
WRITE(6,8040) ISCAL,PMEAN(ISCAL) ,SIGMACISCAL, ISCAL) ,QPVAR, 
$ XMEANI (ISCAL) , COVI(ISCAL, ISCAL) 
NLEFT = NLEFT-1 
210 CONTINUE 
END IF 


IF (QDIAG) THEN 
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CALL PRTMAT(COVI ,NSCALE,CM(KSNAM) ,NSCALE,CM(KSNAM) ,NSCALE, 
$ HEDCOV) 
END IF 


CESTEP 7. 


ACCUMULATE TDM AND COV WITHOUT IMPUTATIONS 
EQUATIONS (1.7) AND (1.8) 


aQqaaaanaa Q °e@ 


DO 340 JSCAL#=1,NSCALE 
DO 330 IRANK=1 ,NRANK 
TCMCIRANK, JSCAL) = TDMCIRANK , JSCAL)+ 
$ WGHT*T (IRANK) *XMEANI (JSCAL) 
330 CONTINUE 
XMFAC#XMEANI ( JSCAL) -PMEAN (JSCAL) 


DO 350 ISCAL=JSCAL,NSCALE 
COV(ISCAL, JSCAL) = COV(ISCAL, JSCAL)+ 
$ WGHT*COVI(ISCAL, JSCAL) 
COVCISCAL, JSCAL) = COV(ISCAL, JSCAL)+ 
$ WGHT* (XMEANI (ISCAL) -PMEAN (ISCAL) ) *XMFAC 
COV(JSCAL ,, ISCAL)=COV(ISCAL, JSCAL) 
350 CONTINUE 
340 CONTINUE 


### END OF CODE TO COMPUTE TDM AND COV 


CESTEP 8. 


### CODE TO COMPUTE THE APPROXIMATE (WEIGHTED) LOG LIKELIHOOD 
COVI IS REUSED TO HOLD A VARIANCE MATRIX 


CESTEP 8(A) 


USE NCAM4AL APPROXIMATION VALUES IF MODE IS NOT FOUND 


IF( CICONV .EQ. 0) )THEN 
QFORM=0 . ODO 
DO 382 JSCAL=1 ,NSCALE 
QFORM=QFORM+( XMEANI(JSCAL)-PMEAN(JSCAL) )*( XMEANI(JSCAL)- 
$ PMEAN(JSCAL) )*SINV(JSCAL, JSCAL) 
QJDIFF=2.0*( XMEANI(JSCAL)-PMEAN(JSCAL) ) 
DO 384 ISCAL=(JSCAL+1) ,NSCALE 
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~3 
fo) 


QFORM=QFORM+QJDIFF*( XMEANI(ISCAL)- 
$ PMEAN(ISCAL) )*SINV(ISCAL, JSCAL) 
? 384 CONTINUE 
382 CONTINUE 
LKTOT=LKTOT+ WGHT*( (XLNLIK-LKNORM)+0.5D0*( DLOG(DETG)-DLSIG ) 
‘ $ -0.SD0*QFORM + D2NSCA*DLOG(2.0D0*PI) ) 


a 


ELSE 
CESTEP 8(B) 


APPLY CORRECTION IN (?7) ACCURATE UP TO N°“(-1/2) ONLY 


agaaan-AaAA Q ° 


QFORM=0.0DO 
DO 392 JSCAL=1,NSCALE 
QFORM=QFORM+( XMODEI(JSCAL)-PMEAN(JSCAL) )*( XMODEI(JSCAL)- 
$ PMEAN(JSCAL) )*SINV(JSCAL, JSCAL) 
QJIDIFF=2.0*( XMODEI(JSCAL)-PMEAN(JSCAL) ) 
DO 394 ISCAL=(JSCAL+1) ,NSCALE 
QFORM=QFORM+QJDIFF*( XMODEI(ISCAL) - 
$ PMEAN(ISCAL) )*SINV(ISCAL, JSCAL) 
394 CONTINUE 
392 CONTINUE 
LKTOT#LKTOT+ WGHT*( (XLNLIK-LKNORM)+0.5D0*( DLOG(DETG)-DLSIG ) 
$ -0.5D0*QFORM + D2NSCA*DLOG(2.0D0*PI) ) 
END IF 


CESTEP 8. END 


### REPLACE GROUP MEAN FROM IMPUTATION ESTIMATE WITH 
GROUP. MEAN FROM ANALYTIC CALCULATION 


aaaaa Q @ 


DO 440 K=1,NTOTG 
ILEV = GVEC(K) 
IF (ILEV .GT. 0) THEN 
DO 430 ISCAL#=1,NSCALE 
GMUCISCAL, ILEV) = GMUCISCAL, ILEV)+WGHT*XMEANI (ISCAL) 

430 CONTINUE 

END IF 
440 CONTINUE 


GO TO 100 
500 CONTINUE 


RETURN -2LOGLIK AND NORM WEIGHTS TO ADD TO SAMPLE SIZE 
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LKTOT=-2D0*LKTOT* (DFLOAT(ISUBJ) /WGTSUM) 
AVEIT=AVEIT/DFLOAT(ISUBJ) 


DO 530 K=1,TOTLEV 
WT = GWN(K) 
IF (WI .NE. ODO) THEN 
DO 520 ISCAL#=1 ,NSCALE 
GMUCISCAL,K) = GMUCISCAL,K)/WT 

520 CONTINUE 

END IF » 
530 CONTINUE 


8000 FORMAT(’ <DIAG> PRINTOUT FOR SUBJECT NUMBER ’,14, 
$ . CYCLE = ’,14) 
8010 FORMAT(’OSCALE’ ,I3,’ THETA LK PRIOR POST’, 
$ , N.POST’) 
8020 FORMAT(’ ’,I5,F10.3,4F10.5) 
8030 FORMAT(’?0 0 ==== = PRIOR-s--99--- eee POSTERIOR----- : 
$ /’ SCALE MU SIGMA VAR(LK) MEAN VAR ’) 
8040 FORMAT(’ ’,15,5F10.5) 
Cc 
Ca## NT FIXING ALTERNATE RETURN CALLS TO SYMINV 
Cc 
7000 WRITE(6,7010) 
7010 FORMAT(’ CALL TO SYMINV FAILED IN ESTEP’) 
STOP 
END 
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Figure 1: Examinee posterior variances for the measurement scale computed using the 
non-modal normal approximation versus the examinee posterior variances computed using 
numerical quadrature with the same values of the parameters input to both procedures. 
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Figure 2: Examinee posterior variances for the measurement scale computed using the 


asymptotic corrections versus the examinee posterior variances computed using numerical 
quadrature with the same values of the parameters input to both procedures. 
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